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Probing Strong Field Gravity Through Numerical 

Simulations 


This chapter describes what has been learned about the dynamical, strong 
field regime of general relativity via numerical methods. There is no rigorous 
way to identify this regime, in particular since notions of energies, velocities, 
length and timescales are observer dependent at best, and at worst are not 
well-defined locally or even globally. Loosely speaking, however, dynamical 
strong held phenomena exhibit the following properties: there is at least one 
region of spacetime of characteristic size R containing energy E where the 
compactness 2GE/c^R approaches unity, local velocities approach the speed 
of light c, and luminosities (of gravitational or matter helds) can approach 
the Planck luminosity jG. A less physical characterization, though one 
better suited to classifying solutions, are spacetimes where even in “well- 
adapted” coordinates the non-linearities of the held equations are strongly 
manifest. In many of the cases where these conditions are met, numerical 
methods are the only option available to solve the Einstein held equations, 
and such scenarios are the subject of this chapter. 

Mirroring trends in the growth and efficacy of computation, numerical 
solutions have had greatest impact on the held in the decades following the 
1987 volume [1] celebrating the 300^^ anniversary of Newton’s Principia. 
However, several pioneering studies laying the foundation for subsequent 
advances were undertaken before this, and they are briehy reviewed in sec¬ 
tion |7.1| below. Though this review focuses on the physics that has been 
gleaned from computational solutions, there are some unique challenges in 
numerical evolution of the Einstein equations; these as well as the basic 
computational strategies that are currently dominant in numerical relativ¬ 
ity are discussed in section 7^ As important as computational science has 
become in uncovering details of solutions too complex to model analyti¬ 
cally, it is a rare moment when qualitatively new physics is uncovered. The 
standout example in general relativity is the discovery of critical phenom- 
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ena in gravitation collapse (section 7.3.1); another noteworthy example is 
the formation of so-called spikes in the approach to cosmological singularities 
(section 7.3.7). A significant motivation for obtaining solutions in the dy¬ 
namical strong field has been to support the upcoming field of gravitational 
wave astronomy, which requires predictions of emitted waveforms for opti¬ 
mal detection and parameter extraction. The expected primary sources are 
compact object mergers, where numerical methods are crucial in the model¬ 
ing of the hnal stages of the events. Binary black hole systems are discussed 
in section 7.3.2, black hole-neutron star and binary neutron stars systems 
in section 7.3.3 Though not of astrophysical or experimental relevance— 
barring the existence of an unexpectedly small Planck energy scale—the 
ultra-relativistic limit of the two body problem is of considerable theoretical 
interest, and this is discussed in section 7.3.5 Spurred by the gauge-gravity 
dualities of string theory, the study of higher dimensional gravity has been 
very active in the past decade; related numerical discoveries are presented in 
section 7.3.6 Some miscellaneous topics are mentioned in section [7. 3. 8[ and 
we conclude the review in section 7.4 with a discussion of open problems for 
the coming years. 

Regarding notation, for the most part we report results in geometric units 
where Newton’s constant G and the speed of light c are set to unity, though 
for clarity some expressions will explicitly include these constants. In re¬ 
ferring to the dimensionality of a manifold, metric or tensor field, we will 
use lower case “d” for spacetime dimensions, and upper case “D” for purely 
spacelike dimensions; e.g., “4d” refers to 3 -|- 1 spacetime dimensions (this 
latter “n 1” form we will also use), and “3D” means three spatial dimen¬ 
sions. 


7.1 Historical Perspective 

The similarly oriented book released in 1987 [I] gave a snapshot of the 
various interesting subjects and problems in gravitational research. However, 
there was no chapter on numerical solution of the Einstein equations, even 
though the subfield of numerical relativity had been in active development 
for over a decade by then. The discipline was still coming into its own, 
and the breadth and scope of works within its purview was still limited. 
Nevertheless, these incipient studies did provide a hint of developments to 
come as the know-how, computational resources and experience improved. 
It is thus important to set some perspective by describing a subset of works 
leading to the current status of the field. 

The particular topics we review later in this chapter are weighted towards 
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developments that have occurred within the past decade or two. This is 
natural as numerical relativity has been a rapidly growing field during this 
time. However, as mentioned, the foundations for building a mature field 
were initiated before, and here we briefly discuss, loosely organized by sub¬ 
ject, some of these more important early results. Unfortunately due to space 
limitations we cannot mention all the relevant works, nor discuss those we 
do mention in any detail. Also, we do not include results, in particular the 
more recent ones, that are discussed elsewhere in this review. 

Binary black hole mergers. The first attempt at a numerical solution 
of the binary black hole merger problem was made by Hahn and Lindquist [2] 
(1964). At that time the term “black hole” had not yet been coined, and the 
full significance of the problem, in particular with regards to gravitational 
wave emission and black hole mergers in the universe, was not recognized. 
Using Gaussian normal coordinates, Misner’s “wormhole” initial data [3], 
representing two black holes initially at rest, was evolved until t k, m/2 
(with m = sJ/AJiQ vr, A being the area of each throat). At that point nu¬ 
merical errors had grown too large to warrant further evolution but it was 
nonetheless possible to measure the mutual attraction between the holes, 
and the fact that the throats were beginning to pinch off. Smarr and 

Eppley m independently revisited the head-on collision problem a decade 
later, now with a profound new understanding of black holes gained in the 
preceding years, both from theory, and from observations suggesting that 
they likely exist in the universe. These works used the same initial data 
as Hahn and Lindquist, Cadez coordinates to simultaneously conform to 
the throats and approach the usual spherical polar coordinates at large dis¬ 
tances [7], and maximal slicing . The culmination of these studies showed 
that the collision emitted radiation of order 0.1% of the total mass, and that 
the waveform was very similar to that computed from a perturbative calcu¬ 
lation (the first indications of the “relative simplicity” of black hole merger 
waveforms discussed in Sec. 


7.3.2). 


In anticipation of construction of the LIGO gravitational wave detectors, 
and the recognized need for waveform models to enable detection, the head- 
on collision calculations were reinitiated by the NGSA group in the early 
’90s [8]. The new simulations offered improved treatment of the Cadez coor¬ 
dinate singularity and radiation extraction, but essentially conhrmed previ¬ 
ous results. With hindsight, it is amusing that in [5] the status of this field 
was summarized as “The two black hole collision problem has been largely 
completed.” This was the prevailing opinion through the mid-’90s, with the 
consensus being that the most significant impediment to solving the full 
3D merger problem was simply lack of available computational power. This 
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turned out not to be the case, and a tremendous effort by the community 
was expended in going from the first short-lived grazing collision simulations 
reported in 1997 [U] to the breakthroughs in 2005 [13 [III [12] that facilitated 
stable evolution of the full problem and the impressive results that have fol¬ 
lowed (see m for more discussion of this development). As briefly discussed 
in Sec. |7.2[ some of the key stumbling blocks were related to the underlying 
mathematical character of the Einstein field equations and the existence of 
geometric singularities inside black holes. This is not to say that limited 
computational power was not an issue; in fact it did hamper the effort to 
rapidly find solutions to the more fundamental problems, as numerous at¬ 
tempts to isolate and solve issues in a symmetry-reduced (or similar) setting 
that could be tackled more quickly with available computational resources, 
failed when carried over to the full problem. 

Gravitational collapse. Numerical studies of the gravitational collapse 
of stars began with the work of May and White [Tl|, who looked at the 
collapse of ideal fluid spheres with a T-law equation of state (specifically 
r = 5/3). They found, depending on the initial conditions, that collapse 
would continue to black hole formation, or halt and then bounce (a neces¬ 
sary condition for an eventual supernova). The “second generation” of codes 
were developed over the next couple of decades, with the pioneering efforts 
of, among others, Wilson US [13, Shapiro and Teukolsky m, Stark and 
Piran |18j . Nakamura |19t I2n| and Evans [2T]. Advances included evolution 
of axisymmetric models to study the effects of rotation and asymmetries, 
solution of the the hydrodynamic equations written in conservative form, 
improvement in the handling of axis coordinate singularities, development 
of moving mesh methods, incorporation of effects of neutrino emission, and 
exploration of a variety of slicing and spatial coordinate conditions. The 


more recent studies of stellar collapse are reviewed in Sec. 7.3.4 


Although many studies of gravitational collapse are motivated by applica¬ 
tion to the wide variety of observed phenomena attributed to stellar collapse, 
there has been considerable work on more theoretical scenarios, in partic¬ 
ular critical collapse, reviewed in Sec. 7.3.1 A notable work we mention 


here is the evolution of collapsing, axisymmetric configurations of collision¬ 
less matter by Shapiro and Teukolsky |22j . In all cases the formation of a 
geometric singularity was observed, but, intriguingly, for sufficiently prolate 
distributions no apparent horizon was found when the time-slice ran into the 
singularity. This could be a slicing issue in that a horizon could still form at 
a later time; however the threshold of prolateness above which no horizons 
were found is consistent with the hoop conjecture [23|, suggesting these cases 
are examples of cosmic censorship conjecture violation in asymptotically flat 
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spacetimes. Another work of theoretical interest that had signihcant impact 
on the foundations of the field was the numerical study of cylindrical gravi¬ 
tational wave spacetimes by Piran [24j (here “collapse” of non-linear waves 
always leads to naked singularities, though the spacetime is not asymptoti¬ 
cally flat). In particular, the modern notions of free vs constrained evolution 
were introduced, and the utility of using coordinate conditions to modify 
the structure of the numerical scheme was demonstrated. 

Binary neutron star, black hole/neutron star mergers. Unlike 
the binary black hole problem which featured extensive early development 
around the head-on collision case, relatively little work on the full general 
relativistic modeling of binary neutron star or black hole/neutron star merg¬ 
ers, head-on or otherwise, was undertaken until the early 2000’s, as reviewed 


in Sec. |7.3.3[ A notable exception is the head-on collision study done by Wil¬ 
son in the late 1970s [16]; he found (by applying the quadrupole formula to 
the matter dynamics) that, similar to the black hole case, ~ 0.1% of the rest 
mass of the spacetime is emitted in gravitational waves. 

Initial Data. The numerical initial data problem in general relativity de¬ 
serves an entire chapter by itself, and unfortunately we are unable to devote 
space to it in this article (for reviews see [25] |26j EZ] , and the books men¬ 
tioned below). We would however be remiss not to mention the York formal¬ 
ism for the construction of initial data |281 |29j . This has become the standard 
method for producing generic initial data for a wide range of problems. More¬ 
over, it provides the framework in which modern ADM-based m Cauchy 
evolution schemes are written, in particular being the starting point to de¬ 
velop the now commonly employed BSSN formalism discussed in Sec. 7.2 A 


few other notable initial data-related works include the Bowen-York closed- 
form solutions to the momentum constraints for black hole initial data ISD, 
the “puncture” initial data |32| (which has had signihcant inhuence beyond 
initial data, leading to the stable evolution of black hole spacetimes without 
the need for excision) and the use of apparent horizons to provide boundary 
conditions and implement singularity excision |33j . 

Miscellaneous. We conclude this historical review by listing a few other 
developments of import to the growth of the held in the 80’s and 90’s. 


Cosmology. We discuss work related to cosmological singularities in Sec. 


7.3.7, and numerical studies of cosmic bubble collisions and local inhomo¬ 


geneities in cosmology in Sec. 7.3.8, but mention here that much pioneer¬ 
ing work on numerical cosmologies and spacetimes with related symme¬ 
tries began with the works of Centrella, Anninos, Wilson, Kurki-Suonio, 
Laguna and Matzner, Berger and Moncrief [Ml EHESIETI EH] 
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Boson stars. The numerical study of these self-gravitating soliton-like con- 
hgurations of scalar helds has a long history that begins with calcula¬ 
tions in the late 1960’s by Kaup [32] and by Ruffini and Bonnazola m, 
who found spherically-symmetric, static solutions in the Einstein-Klein- 
Gordon model. A major resurgence of interest in the subject was sparked 
by Colpi, Shapiro and Wassermann’s discovery that the addition of a 
non-linear self-interaction could lead to boson star masses that, in con¬ 
trast to those originally constructed, were in an astrophysically interesting 
range m- Much subsequent work investigating a wide variety of types of 
boson stars and related objects has been carried out since, and we touch 
on some representative calculations in Secs. 7.3.1 and 7.3.5 Once more, 
however, space limitations preclude a thorough coverage of this topic and 
we direct the interested reader to reviews such as 


• Excision. The hrst successful simulation incorporating the use of black 
hole excision to eliminate geometric singularities from the computational 
domain was presented by Seidel and Suen [i3] . 


Hyperbolic evolution sehemes. One of the influential efforts predating the 
wave of activity searching for stable hyperbolic evolution schemes dis¬ 
cussed in Sec. 7.2 was the formulation of Bona and Masso 


• The Grand Challenge (1993-1998). This large scale NSF-funded project, 
aimed at solving the black hole inspiral and merger problem, involved es¬ 
sentially all US-based numerical relativists and, crucially, many computer 
scientists. The notable results culminating from this effort include propa¬ 
gation of a single Schwarzschild black hole through a 3D mesh |l5|, early 
efforts in refined gravitational wave extraction methods and improved 
outer boundary conditions [36], and the development of a characteris¬ 
tic code that could stably evolve even highly perturbed single black hole 
spacetimes m 


7.2 Numerical Relativity: Current State of the Art 

Before beginning our review of the important physics garnered from numer¬ 
ical solutions of the Einstein field equations over the past few decades, we 
describe some of the key insights obtained along both formal and numerical 
fronts that have made these ventures possible. Of course, it is impossible 
to exhaustively cover all of them; we thus choose particularly relevant ones 
that have had a strong influence on the field. 
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7. 2.1 Mathematical Formalism 

Any numerical study of a dynamical process requires solving a suitably for¬ 
mulated initial value (or initial boundary-value) problem. That is, provided 
a set of evolution equations, together with a given state of the system at an 
initial time, its future evolution can be obtained by a numerical integration. 
At face value any covariant theory is at odds with this requirement, unless 
a suitable ‘time’ foliation of the spacetime is introduced. In the case of the 
Einstein equations, projections (tangential and normal to each leaf of the 
foliation) provide a natural hierarchy of evolution and constraint equations. 
The latter are tied to the fact that the Einstein equations are overdetermined 
with respect to the physical degrees of freedom, and allow one to chose dif¬ 
ferent combinations of equations to solve for the spacetime to the future 
of the initial hypersurface. As a result, one can distinguish free evolution 
approaches—where only evolution equations are employed to this effect— 
from constrained (partially constrained) approaches where (some of the) con¬ 
straints are used to solve for a subset of variables [Mj • We also note that the 
freedom in choosing a foliation gives rise to Cauchy (or D-|-l), characteristic 
or hyperboloidal formulations, in the case of spacelike, null or spacelike-but- 
asymptotically-null foliations respectively. For a review of related numerical 
approaches, see e.g., [H]; a more pedagogical exposition of the basic con¬ 
cepts can be found in recent textbooks on the subject ggEolEIlES]. 

On the Cauchy front, early efforts employed the so-called York-ADM for¬ 
mulation [25] (a reformulation of the standard Hamiltonian-based ADM ap¬ 
proach). This formulation is geometrically appealing in that it provides evo¬ 
lution equations for the intrinsic and extrinsic curvatures of the foliation. 
However, beyond spacetimes where symmetries allow for reducing the dimen¬ 
sionality of the problem, numerical evolution with the York-ADM method 
exhibits instabilities. In the early 2000s [53], it was recognized that such a 
formulation is only weakly hyperbolic, implying that at the analytical level 
the system of PDFs lack properties required to achieve robust numerical im¬ 
plementation [5l]. A flurry of activity in the following years provided much 
insight on how to deal with this issue and construct (desirable) symmetric 
or strongly hyperbolic formulations by suitable modifications of the equa¬ 
tions (primarily via the addition of constraints and the use of appropriate 
coordinate conditions—see e.g. [55] [56] [57] for reviews on the subject). It 
turns out that arbitrarily many formulations could be dehned with these 
desirable properties and all are, of course, equivalent at the analytical level. 

At the discrete level however, this is not the case. In particular in free evo¬ 
lution schemes it is challenging to control the magnitude of the truncation 
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error associated with the constraint equations (which are not explicitly im¬ 
posed except at the initial time). Such errors compound at different rates in 
different formulations, and the practical physical time that accurate results 
can be achieved in simulations thus varies significantly. Two formulations 
have been shown empirically to display the robustness needed to construct 
a large class of 4-dimensional spacetimes (without symmetries) that are 
of particular relevance to the contemporary astrophysical and theoretical 
physics problems discnssed here. These are the generalized harmonic evolu¬ 
tion with constraint damping [10] (closely related to the Z4 formalism, see 
e.g. |58|), and the BSSN (or BSSNOK) approach [5^[60ll6T] . 

Harmonic coordinates have a history even older than the field equations 
themselves—having been used by Einstein in his search for a relativistic 
theory of gravity as early as 1912 |62|— and have played an important role 
in many key discoveries of properties of the field equations since (see the in¬ 
troduction in |63j ). Harmonic coordinates can be defined as the requirement 
that each spacetime coordinate obeys a homogeneons scalar wave equation. 
Enforcing this at the level of the Einstein equations converts the latter to 
a form that is manifestly symmetric hyperbolic. This desirable property is 
maintained if freely specifiable functions are added as source terms in the 
wave equations, resulting in generalized harmonic coordinates. In principle, 
the source functions allow arbitrary gauges to be implemented within a 
harmonic evolution scheme [M]. Constraint damping terms m are further 
added to tame what otherwise would result in exponential growth of trun¬ 
cation erroi[^ This formulation has also been particnlarly nseful in finding 
stationary black hole solutions in higher dimensions where, due to the sta¬ 
tionary nature of the spacetime, the coordinate freedom can be exploited 
to define a convenient, strictly elliptic problem |68j . The BSSN formulation 
is an extension of the York-ADM approach which introdnces several addi¬ 
tional (constrained) variables to remove particnlar offending terms from the 
equations, and this, together with a judicious choice of coordinates, ensures 
strong hyperbolicity of the underlying equations. 

When black holes are evolved, the geometric singnlarities inside the hori¬ 
zons need to be dealt with in some manner to avoid numerical problems 
with infinities. One such method is excision |33j . where an inner excision 
boundary is placed inside each apparent horizon to remove the singular re¬ 
gion from the computational domain. Due to the causal structure of the 

^ The idea of constraint damping via the addition of terms to the equations that are 

homogeneous in the constrained variables (hence are zero for continuum solutions) dates back 
several years earlier for the Einstein equations m, and since has also effectively been applied 
to other systems of PDEs that have internal constraints, in particular Maxwell’s 
equations EH. 
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spacetime, the characteristics of the evolution equations ostensibly all point 
out of the computational domain at the excision boundary, and no bound¬ 
ary conditions are placed there. Excision is commonly used in harmonic 
evolution schemes. For evolution within the BSSN approach, an alternative 
moving puncture method has proven successful mm- This is an extension 
of puncture initial data, where the puncture point inside the horizon that 
formally represents spatial infinity on the other side of a “wormhole” is now 
evolved in time. With the typical gauges employed during evolution the in¬ 
terior geometry evolves to a so-called “trumpet” slice, where the puncture 
asymptotes to the future timelike infinity of the other universe [69l [70] . Effec¬ 
tively then, the puncture also excises the singularity from the computational 
domain. 

On the characteristic front, the structure of the equations is significantly 
different from the Cauchy problem, as the foliation is defined using char¬ 
acteristic surfaces. The system of equations displays a natural hierarchy of 
evolution equations, constraints, and a set of hypersurface equations for vari¬ 
ables that asymptotically are intimately connected to the physical degrees 
of freedom of the theory m- Beyond spherically symmetric applications, 
numerical codes employing this formulation show a remarkable degree of 
robustness, and can stably evolve highly distorted single black hole space- 
times mm- However, the rigidity in the choice of coordinates (being tied 
to characteristics) imply that difficulties arise when caustics and crossovers 
develop. For astrophysical purposes, the main role of the characteristic ap¬ 
proach has been to provide a clean gravitational wave extraction proce¬ 
dure mm and to study isolated black holes. Outside of the astrophysical 
domain, it has been convenient in studies of black hole interiors (e.g. m) 
and has become the predominant approach used to exploit the AdS/CFT 
(Anti-deSitter/Conformal Field Theory) duality of string theory (see [TH] 
and references therein). Here gravity in asymptotically AdS spacetimes is 
used to study field theory problems, some examples of which are described 
in section 17.3.61 

Finally, the hyperboloidal formulation m adopts a Cauchy approach in 
a conformally related spacetime where the physical spacetime is recovered 
as a subset of a larger one. This allows studying, within a single framework, 
both the local and global structure of spacetime (as the larger manifold 
covers spacelike, null and timelike infinity in a natural way). However, it has 
received considerably less attention than the other approaches, though some 
interesting first steps have been carried out (e.g. [781I79]). 
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7.2.2 Numerical Methods 

The subject of numerical analysis as it pertains to solutions of problems in 
applied mathematics is, of course, vast. Even when restricting to what is 
relevant for numerical relativity applications, the breadth of methods em¬ 
ployed is considerable and, naturally, depends on the particular goal one 
has in mind. From constructing initial data and evolving the solution to the 
future of an initial hypersurface, to extracting physical information from the 
numerical results, an abundance of different techniques and methods have 
been used. Here we provide some brief comments with an aim to impart a 
basic understanding of the available options. 

Any numerical implementation ultimately renders the problem of interest 
into an algebraic problem for a discrete number of variables that describe 
the sought-after solution. For gravitational studies, this involves devising ap¬ 
proximate methods to solve the relevant partial differential equations. Such 
methods can be conveniently visualized as providing a way to discretize the 
underlying variables that describe the problem as well as the spatial deriva¬ 
tives within the equations, and providing a recipe to advance the solution 
in discrete time. 

The technique most commonly used is the finite difference (FD) method, 
with the solution represented by its value at discrete grid-points covering 
the manifold of interest. Discrete spatial derivatives are defined through 
suitable Taylor expansions, which can allow for high order accurate approx¬ 
imations for smooth solutions. Further refinements can be achieved through 
the use of discrete derivative approximations that satisfy summation by 
parts. This property is a direct analog of integration by parts often exploited 
to obtain estimates about the behavior of general solutions at the analyt¬ 
ical level [52 [55]. What has proven especially useful for many problems 
is the adoption of adaptive mesh refinement (AMR) adaptive mesh refine¬ 
ment to efficiently resolve a large range of relevant spatio-temporal length- 
scales, and without a priori knowledge of the development of small scale 
features ini EH E2|. Another discretization approach used is the pseudo- 
spectral method [83l [82 : here the solution is expanded in terms of a suitably 
chosen basis (e.g. Chebyshev or Fourier), which also provides a simple way 
to compute spatial derivatives. The coefficients of the expansion provide the 
sought after solution. Pseudo-spectral methods provide a highly efficient 
way to achieve high accuracy results for smooth solutions. As with AMR for 
FD, adaptive, multi-domain decomposition methods can be employed for 
efficiently resolving the length scales in the problem (see e.g. |85jl. 

Advancing the solution in time requires integrating the discrete values of 
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the FD solution, or spectral coefficients, via suitable operators. One common 
approach is the method of lines, where having discretized the spatial part 
of the problem, accurate methods devised for ordinary differential equations 
are used to integrate the variables in time. 

The above discretization approaches are also well-suited to evolving ad¬ 
ditional fields coupled to gravity that are smooth and, in particular, do not 
develop discontinuities (such as scalar or electromagnetic fields). However, 
when matter sources such as neutron stars are incorporated, the equations of 
general-relativistic hydrodynamics (or magnetohydrodynamics) must also be 
solved. These can be expressed in a way fully consistent with the approaches 
employed to integrate Einstein equations (for a review on this topic see [86] ). 
Nevertheless, as the solution to the hydrodynamic equations can induce dis¬ 
continuities even for smooth initial conditions, finite volume methods m 
are most often adopted as they are especially suited to accurately handle 
such features. 


7.3 Strong Field Gravity 

As mentioned, numerical simulations are often the only way to gain insights 
into the behavior of gravity in the strongly non-linear regime. This arises 
naturally in gravitational collapse, systems involving black holes and neu¬ 
tron stars, ultra-relativistic collisions and cosmology. Here questions range 
from fundamental explorations of the theory itself, to the resolution of ques¬ 
tions of astrophysical relevance, such as the characteristics of gravitational 
wave signals produced in compact binary mergers, or the effect of highly 
dynamical and strong gravitational fields on matter/gas/plasma and their 
role in powering spectacular phenomena like gamma ray bursts. Complex 
numerical simulations have been developed over the past decades to start 
answering these long-standing questions, and have produced results that 
often raise new questions. In what follows, we discuss some of the more im¬ 
portant findings and open questions, organizing the presentation by subject 
area. These examples are necessarily limited in scope and presentational 
depth, but serve as illustrations of the breadth of problems addressed with 
simulations. For a complimentary review article on numerical relativity and 
its applications, see [88] . 
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7.3.1 Critical Phenomena in Gravitational Collapse 

Overview 

Gravitational collapse, including the process of black hole formation, is one 
of the hallmarks of general relativity. As has already been noted in Sec. |7.1[ 
although simulations play an increasingly dominant role in advancing our 
understanding of collapse scenarios that we believe play out in the universe, 
they also provide the means to perform detailed studies of more fundamen¬ 
tal aspects of the process. Albeit reflective of more than a little theorist’s 
conceit, we can view computer programs as numerical laboratories which— 
paralleling real experiments in nonlinear science—are endowed with one or 
more control parameters that are varied in order to unearth and elucidate 
the phenomenology exhibited by the setup. 

Critical phenomena are concerned with families of solutions to the coupled 
Einstein-matter equations (including the vacuum case), where a continuous 
parameter p labels the family members, and it is assumed that the space- 
times are constructed dynamically—usually via simulation—starting from 
prescribed initial data that depends on p. The initial data typically rep¬ 
resents some bounded distribution of initially imploding matter and p is 
chosen to control the maximal strength of the gravitational interaction that 
ensues. For p sufficiently small, gravity remains weak during the evolution, 
and the spacetime is regular everywhere (if the matter is massless radiation, 
for example, the radiation will disperse to infinity, leaving flat spacetime in 
its wake). For p sufficiently large, gravity becomes strong enough to trap 
some of the matter in a black hole, with mass MbH) and within which a sin¬ 
gularity forms. For some critical value p* lying between the very-weak and 
very-strong limits, the solution corresponds to the threshold of black hole 
formation, and is known as a critical spacetime for the given model. Collec¬ 
tively, the properties of these special configurations, as well as the features 
associated with the spacetimes close in solution space to the precisely criti¬ 
cal solution, comprise what is meant by critical behaviour. Evidence to date 
suggests that virtually any collapse model that admits black hole formation 
will contain critical solutions. 

It transpires that Mbh can be formally viewed as an order parameter, in 
the statistical mechanical sense, and most of the critical solutions identified 
to date can be sorted into two basic classes based on the behaviour of Mbh 
at threshold. Specifically, solution-space behaviour corresponding to both 
first- and second-order phase transitions is seen, defining what are called 
Type I and Type II critical solutions, respectively. Thus, in the Type I case 
A/bh is finite at the threshold, so Mbh(p) exhibits a gap (jump) at p = p*. 
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Conversely, in the Type II instance Mbh becomes infinitesimal as p —)■ p* 
from above, and there is no gap. It is also crucial to observe that the precisely 
critical solution for p = p* does not contain a black hole. 

There are three key features associated with both types of black hole 
critical phenomena: universality, symmetries and scaling. 

Concerning the first property, most, if not all. Type II solutions, and 
some Type I solutions, exhibit a type of universality in the sense that one 
finds the same critical configuration through numerical experimentation as 
sketched above, irrespective of the specific way p parameterizes the initial 
data. This implies a certain type of uniqueness, or at least isolation, of the 
critical spacetimes in solution space, analogous to the uniqueness of the 
Schwarzschild solution as the endpoint of black hole formation in spherical 
symmetry. 

Secondly, Type I critical solutions generically possess a time-translational 
symmetry, which is either continuous (so the solution is static or stationary), 
or discrete (so the solution is periodic). In the discrete case, the oscillation 
frequency forms part of the precise description of the critical solution, and 
is usually determined by an eigenvalue problem. Type II critical solutions, 
on the other hand, typically have a scale, or homothetic, symmetry and are 
therefore scale invariant. Once again, this symmetry can be either continuous 
or discrete (CSS/DSS for continuously/discretely self-similar). Continuously 
self-similar solutions have long been studied in relativity as well as in many 
other areas of science, frequently arising in situations where the underlying 
physics has no intrinsic length scale. On the other hand, it is safe to say 
that the observation of discrete self-similarity in the earliest numerical cal¬ 
culations of critical collapse came as a complete surprise. For DSS solutions, 
the analogue of the frequency of periodic Type I configurations is known 
as the echoing exponent, A. When expressed in coordinates adapted to the 
self-similarity, a DSS solution is oscillatory with period A; each complete 
oscillation represents a shrinking of the scale of the dynamics by a factor of 
e^. In relation to the collapse process, self-similar behaviour of either type 
is particularly interesting because it means that at criticality the strong field 
regime propagates to arbitrarily small spatiotemporal scales. Indeed, as the 
self-similar solution “focuses in” to an accumulation event at the center 
of the collapse, curvature quantities grow without bound, but with no for¬ 
mation of an event horizon. Thus, Type II critical solutions possess naked 
singularities and have significant relevance to the issue of cosmic censor¬ 
ship. However, in terms of their origin from collapse, it is imperative to note 
that these naked singularities are produced from (infinite) fine tuning of the 
initial data, so are therefore not generic with respect to initial conditions. 
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The third characteristic of critical collapse is the emergence of scaling 
laws. Empirically, these are determined by studying the behaviour of certain 
physical quantities as a function of the family parameter p as p —)• p*. For 
Type I behaviour a typical scaling law measures the time interval, r, during 
which the dynamically-evolving conhguration is close to the precisely critical 
solution. One finds 

T ~—cr In Ip — p*| (7.3.1) 

where a is called the time-scaling exponent. In Type II collapse the black 
hole mass scales according to 

MBH~C|p-p*r (7.3.2) 


where 7 is known as the mass-scaling exponent, and C is a family dependent 
constant. For both types of behaviour, if the critical solution is universal with 
respect to the initial data, then so is the corresponding scaling exponent. 
Again, this is the case for all known Type II solutions, but not so for most 
Type I transitions where, as discussed in more detail below, a particular 
critical configuration is typically one member of an entire branch of unstable 
solutions. 

The scaling laws can be understood in terms of perturbation theory. 
The key observation [89l [901 E] is that the appearance of critical solu¬ 
tions through a fine tuning process—wherein one of two distinct end states 
characterizes the long-time dynamics—suggests that they have a single un¬ 
stable perturbative mode. The inverse of the Lyapunov exponent associated 
with that mode is then precisely the scaling exponent. Furthermore, leading 
subdominant modes can give arise to additional scaling laws, for example 
charge and angular momentum 


The scaling relations (7.3.1)-(7.3.2) underscore the fact that near criti¬ 
cality, there is exponentially sensitive dependence on initial conditions, and 
that irrespective of the original choice of initial data parametrization, it is 
the transformed quantity In \p — p*\ which is most natural in describing the 
phenomenology. In simulations one wants to compute with \p — p*\ as small 
as possible in order to most accurately determine the threshold solutions and 
their associated exponents. This is especially true for the Type II case since 
the asymptotically flat boundary conditions that are normally adopted are 
incompatible with self-similarity, so one relies on calculations which probe 
as small scales as possible to ensure that boundary effects are minimized. 
In practice, investigation of the critical regime is ultimately limited by the 
fact that p can only be fine-tuned to machine precision. For spherical and 
axisymmetric calculations, this can be accomplished with current technolo- 
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gies, but, for the Type II scenarios, only if the numerical algorithm provides 
sufficient spatio-temporal dynamic range using a technique such as adaptive 
mesh-refinement (AMR) [SU]. Indeed, at this point it should be emphasized 
that almost all known critical solutions have been computed in models that 
impose symmetry restrictions. Spherical symmetry has been most commonly 
adopted. There have been some axisymmetric calculations, but very few fully 
3D studies. 

As already mentioned, the appearance of critical solutions seems com¬ 
pletely generic, irrespective of the matter content of the model. However, 
details of the phenomenology are dependent on a number of factors in¬ 
cluding the following; the type of matter and the specifics of any self¬ 
interaction terms, the imposed symmetries, the spacetime dimensionality 
and the asymptotic boundary conditions. The remainder of this section is 
devoted to a summary of a necessarily incomplete selection of the many nu¬ 
merical studies performed to date, organized by the type of matter employed. 
In order to highlight the state of the art in the subject, there is some bias 
towards more recent calculations, and an attempt has been made to impart 
some sense of the wide variety of scenarios that have been explored. Those 
interested in more information are directed to the excellent comprehensive 
reviews of the subject [Ml E2]- 


Scalar Fields 

Critical collapse was first studied in the model of a spherically symmetric 
minimally coupled massless scalar field [9Tj. Using several families of initial 
data, a single Type II DSS solution with 7 « 0.37 and A « 3.44 was found 


(see Fig. 7.1). Due to the extreme dynamical range required to fully resolve 


the critical solution, use of adaptive mesh refinement was crucial. Remark¬ 


ably, the mass-scaling relation (7.3.2) provided a good fit for Mbh even when 
\p—p*\ was large enough that the final black hole contained most of the total 
mass of the spacetime. The main results of [93] have been confirmed many 


times since using a variety of different algorithms and coordinate systems. 
Assuming the existence of a spacetime with a discrete homotheticity, certain 
regularity conditions and a tailored numerical approach, the critical solution 
and associated exponents were computed to very high accuracy in [95] . Ad¬ 
ditionally, analysis of the implications of the discrete self-similarity led to the 


prediction and observation of a modulation of the mass-scaling law (7.3.2) 
with period A/( 27 ) |95l I96j . Finally, an important analysis in |97] found 
that all non-spherical modes of the critical solution decay, strongly suggest¬ 
ing that the same threshold configuration would appear if the symmetry 
restriction were relaxed. 
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Figure 7.1 Type 11 discretely self-similar critical solution computed from 
the collapse of a spherically symmetric distribution of massless scalar field. 
The figure shows the late time configuration of the scalar field from a 
marginally subcritical evolution where the family parameter has been tuned 
to approximately a part in 10^®. The radial coordinate is logarithmic, mak¬ 
ing the discretely self-similar (echoing) nature of the solution evident: each 
successive echo represents a change in scale of e‘^ « 31. The data were 
generated using the axisymmetric code described in |98j . 


Critical solutions from axisymmetric massless scalar collapse using mul¬ 
tiple initial data families were constructed in [98]. For the most part the 
threshold configurations could be described as the spherical solution plus 
perturbations (measured values for the scaling exponents were 7 ~ 0.28- 
0.41 and A ~ 2.9-3.5), but there were also indications of a single asymmet¬ 
ric mode which, as it grew, produced two separated regions (on axis) within 
which the solution locally resembled the spherical one. This observation is 
in conflict with IHZ!, bnt the accuracy of the results was insufficient to con¬ 
vincingly demonstrate that the growth was gennine and not a reflection of 
limitations in the simulations. Adaptive mesh refinement was again crucial. 

Very recently, a stndy of massless scalar collapse using a fully 3D code has 
been carried out [99| and, in fact, represents the first calcnlations of Type 
II general-relativistic critical phenomena in 4 spacetime dimensions without 
symmetry restrictions. Fonr initial data families defining a spherical mat¬ 
ter distribntion deformed to varying degrees with a Y 21 spherical harmonic 
anisotropy were considered. Even though AMR was used, compute-time lim¬ 
itations kept the tuning of p/p* to about a part in 10^. Nonetheless, evi- 
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dence for the emergence of the spherically symmetric critical solution with 
7 ~ 0.37-0.38 was found. There were also preliminary indications of echoing 
with A ~ 3.1-3.3. 

The massless scalar model has no intrinsic length scale so, in retrospect, 
the appearance of a Type II solution at threshold is natural. Introduction 
of a mass, fi, breaks scale invariance and, as shown in [lOOj . complicates the 
picture of criticality. For initial data with a length scale A the massless be¬ 
haviour is recovered when A/i <C 1. However, for A/x > 1, a Type I transition 
is seen with a critical solution which is one of the periodic, star like configu¬ 
rations (oscillons) admitted by the model and constructed in |101j . As with 
relativistic perfect fluid stars, the oscillons comprise a one-parameter family 
that can be labeled by the central density. As the central density increases 
the stellar mass also increases, but only up to a point, whereafter dynamical 
instability sets in and the stars reside on the so-called unstable branch—it 
is precisely one of these unstable solutions that sits at the Type I transition. 
This latter type of behaviour was also observed in m using a massive 
complex scalar field whose static solutions, known as boson stars, also have 
stable and unstable branches. In this instance stable stars were driven to 
a Type I threshold via an imploding pulse of massless scalar field, whose 
overall amplitude was used as the tuning parameter. 

Investigation of circularly symmetric massless scalar collapse in 2 -|-1 AdS 
spacetime [10311104) represents one of the few instances where critical be¬ 
haviour in a non-asymptotically flat setting has been seen (but also see 
the discussion of the turbulent instability of 3 -|- 1 AdS [105) in Sec. 7.3.8). 
Evidence for a Type II transition with a CSS solution was found—with a 
mass-scaling exponent 7 ss 1 . 2 —but a thorough understanding of the pic¬ 
ture of criticality here is still lacking. In particular, an analytic CSS solution 
that shows good agreement with the numerical results has been found |106j . 
but seems to have additional unstable modes. Its existence also seems para¬ 
doxical in the sense that, heuristically, the cosmological constant should be 
irrelevant on the small scales pertinent to scale-invariance, yet is essential 
in the construction of the solution. 


Vacuum 

Historically, the second example of black hole critical phenomena discovered 
was in the collapse of pure gravitational waves [107) in axisymmetry. The 
study employed one family of initial data representing initially incoming 
pulses of gravitational radiation with quadrupolar angular dependence, and 
with an overall amplitude factor serving as the control parameter. Evidence 
for a Type H transition was found, with a discretely self-similar critical 
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solution that was centred in the collapsing energy. The computations yielded 
an estimated mass-scaling exponent 7 ~ 0.37 and an echoing factor A ~ 0.6. 
The calculations did not use AMR, but due to the use of spherical polar 
coordinates, increased central resolution could be achieved with a moving 
mesh technique. Nonetheless, the dynamic range of the code was very limited 
relative to that used in [9l], so it was quite fortuitous that A in this case 
was quite small. 

It is truly remarkable that in the two decades that have elapsed since the 
publication of |107j . and despite several additional assaults on the problem 
and a vast increase in the available amount of computer resources, little 
progress has been made in reproducing and extending these early results. 
One notable exception is [lOSj in which the collapse of axisymmetric Brill 
waves was studied, using several different families of data with varying de¬ 
grees of anisotropy. Once more, evidence for a Type II transition was found in 
all of the experiments, with a scaling exponent 7 —measured in this instance 
through the scaling of a curvature invariant in subcritical collapse |109j —in 
the range 0.37-0.4. However, in stark contrast to the observations in |107| . 
most of the computed critical solutions showed accumulation on rings at fi¬ 
nite distances from the origin, rather than at the origin itself. Additionally, 
indications of echoing were seen, but with an estimated A ~ 1.1 significantly 
different from that reported in |107) . Development of a more complete un¬ 
derstanding of the critical behaviour of collapsing gravitational waves, both 
in axisymmetry and the full 3D case, remains one of the most important 
unresolved issues in this field. 

In D -|- 1 dimensions, with D even, application of a co-homogeneity two 
symmetry reduction to the vacuum Einstein equations yields a set of wave 
equations dependent only on a single radial dimension. In contrast to those 
resulting from a spherically symmetric reduction, these equations admit 
asymptotically-flat, radiative solutions [llOt 1111] I112t fTT3] . For D = 4, and 
adopting the so-called biaxial ansatz. Type II DSS behaviour was found, 
with A PS 0.47 and 7 ps 0.33 m- Analogous results were found for D = 8 
where A ps 0.78 and 7 ps 1.64 m- The more general triaxial ansatz for 
D = 4 was considered in m- Here, the biaxial critical solution still ap¬ 
pears at threshold, but due to a discrete symmetry in the model, the crit¬ 
ical surface actually contains three copies of the configuration. As well, on 
the boundaries of the basins of attractions of these copies, a different DSS 
solution with two unstable modes was predicted and computed using a two- 
parameter tuning process. Additional numerical experiments have shown 
that the critical-surface boundaries have a fractal structure m • 
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Studies of critical behaviour with perfect fluid sources have been extremely 
important in the development of the subject, not least since it was in this 
context that understanding of the phenomena in terms of unstable pertur¬ 
bative modes was developed. The first calculations focused on spherically 
symmetric simulations with a fluid equation of state (EOS), P = kp, where 
P and p are the fluid pressure and energy density, respectively, and with the 
specific choice /c = 1/3 (radiation fluid) |89j . A continuously self-similar crit¬ 
ical solution was found with a mass-scaling exponent 7 ~ 0.36. In addition, 
the critical solution was computed independently by adopting a self-similar 
ansatz, and was shown to be in excellent agreement with the simulation 
results, and it was suggested that a perturbation analysis could be used to 
at least approximately compute 7 . Such an analysis was carried out in |90|, 
where both the critical solution and its linear perturbations were deter¬ 
mined, and it was shown that there was a single unstable mode whose in¬ 
verse Lyapunov exponent yielded the same value of 7 seen in the simulations. 
Interestingly, at this time the values of 7 that had emerged from the three 
models for which threshold solutions had been identified were numerically 
the same to the estimated level of numerical accuracy, suggesting that the 
mass-scaling exponent might be universal across all matter models. How¬ 
ever, the results of [91] (performed at the same time as [9^), where critical 
solutions and their perturbative modes were determined via the self-similar 
ansatz for many values of the EOS parameter k in the range 0.01-0.888, 
showed definitively that 7 was in general model-dependent. A more exten¬ 
sive perturbation analysis |II4l III5j suggested that the spherical solutions 
will appear at threshold when spherically symmetry is relaxed only for values 
of k in the range 1/9 < /c < 0.49; for other values of k, additional unsta¬ 
ble modes were found. These conclusions have yet to be verified through 
simulations, and it will be very interesting to do so. 

The P = kp EOS is scale-invariant (and in fact is the only EOS com¬ 
patible with self-similarity [US]) so Type II critical behaviour is expected. 
Eor more general equations of state, including the commonly adopted ideal 
gas law, intrinsic length scales appear and, as anticipated from the massive 
scalar field studies |I00j . the critical phenomenology becomes richer. In par¬ 
ticular, the expectation that unstable stars can appear as Type I critical 
solutions was confirmed in |II7j using the ideal gas EOS and the same type 
of experiments performed in [Mj. Type II behaviour with this EOS also 
appears when the fluid internal density of the configuration is much larger 
than the rest energy density |II 8 l III9( 1120] . in which case the EOS limits 
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to the scale-free equation, and the measured mass-scaling exponents, agree 
with those computed from a scale-invariant ansatz. 

A possible cosmological application of Type II fluid collapse was posited 
in m, where it was argued that the mass-scaling relation should apply to 
the formation of primordial black holes, since the exponential decay of the 
scale of density fluctuations entering the horizon at any epoch provides an 
intrinsic fine tuning mechanism. This leads to a modification of the usual 
mass function for the primordial black holes, which incorporates the predic¬ 
tion that holes of sub-horizon scale could form at all times. 

Over the past few years, significant progress has been made in extending 
the investigations of Type I critical behaviour with fluids to the axisymmet- 
ric |122( 1123] I124t fT25] and 3D |126) arenas. Almost all studies have adopted 
a stiff {k = 1) ideal gas EOS (with static or stationary solutions interpreted 
as neutron stars), and the work reported in [126j also incorporated rotational 
and magnetic effects. In |122j a Type I transition was observed in the head- 
on collisions of two neutron stars where several different tuning parameters, 
including the stellar mass and the index k, were employed. Clear evidence 
of lifetime scaling for subcritical evolutions was seen. It was also suggested 
that the change in the EOS that occurs as a real post-collision remnant cools 
could provide a natural tuning mechanism, so that if the cooling was suffi¬ 
ciently slow, the critical solution might have astrophysical relevance. Eurther 
simulations of head-on collisions [12311125| have corroborated these findings, 
and it was demonstrated in |123j that the end state of the marginally sub- 
critical collision was well-described by a perturbed star on the stable branch. 
Intriguingly, the lifetime scaling measured in |123j exhibits a periodic mod¬ 
ulation of cr—analogous to that seen in the mass-scaling exponent for DSS 
Type II transitions—that has yet to be explained. The fact that stars on an 
unstable branch can be identified as Type I solutions has also been demon¬ 
strated in a more direct fashion, through the use of initial data families 
where the tuning parameter perturbs (or effectively perturbs) a star known 
or suspected to be one-mode unstable. This strategy was employed in |124j 
to demonstrate the criticality of an unstable spherical configuration, with 
an accurate computation of cr. Einally, in |126j evidence for the threshold 
nature of rotating unstable stars—both non-magnetized and magnetized— 
with preliminary evidence of lifetime scaling was reported. This last study, 
along with [99j . provides a tantalizing glimpse of what lies in store for this 
field as symmetry restrictions are relaxed and the physical realism of models 
is enhanced. 
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Other Types of Matter 
Spherical collapse of an SU{2) Yang-Mills field within the magnetic ansatz 
was studied in |127| , and was the hrst model where both Type I and Type II 
behaviour was observed. Here the n = 1 member of the Bartnik and McKin¬ 
non countable sequence of static configurations |128j . which had previously 
been shown to have one unstable mode, is the attractor for the Type I tran¬ 
sition, while a DSS solution with A = 0.74 and 7 = 0.20 was also found. 
The model exhibits another transition, strictly in the black-hole sector of 
solution space, where colored black holes arise at the threshold and where 
Mbh has a gap as one tunes across it |129j . 

Spherically symmetric self-gravitating cr-models (wave-maps), which typi¬ 
cally incorporate dimensionless tunable coupling constants, have been shown 
to display especially rich critical phenomenology. Notably, transitions be¬ 
tween CSS and DSS Type II behaviour as the coupling is varied have 
been seen in both the 2-dimensional nonlinear model |130j and the SU{2) 
case m- The transition in the latter instance is particularly interesting, 
displaying behaviour where near-critical evolutions approach and depart 
from a CSS solution episodically. 

Finally, Type I critical behaviour has been seen in the collapse of collision- 
less matter in spherical symmetry—with or without a particle mass—where 
the threshold solutions are static |1321113311134j and appear to exhibit the 
expected properties of Type I solutions, including lifetime scaling. In the 
massless case it has been argued that there should be no one-mode unstable 
solutions m, and this apparent contradiction with the numerical results 
remains another unsolved puzzle. 

7.3.2 Binary Black Hole Mergers 

The non-linear nature of general relativity has several interesting conse¬ 
quences for how it describes particles and the gravitational interaction be¬ 
tween them. First, technical caveats aside, the simplest possible solution 
describing the geometry of an idealized point-like distribution of chargeless, 
spinning matter is a Kerr black hole. From an external observer’s perspec¬ 
tive there is thus no geometrical realization of a point-like structure, as the 
event horizon prevents length scales smaller than the energy (in geometric 
units) of the black hole from being probed. Second, there is no analogue of 
a Newtonian potential that can be superposed to come up with a simple 
description of the interaction of two black holes. In consequence, a detailed 
understanding of one of the most basic interactions in gravity, the two body 
problem, reqnires numerical solution of the field equations. On the other 
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hand, thanks to the “no-hair” properties of black holes, the merger of two 
Kerr black holes is expected to describe the merger of all astrophysical black 
holes essentially exactly, the only idealization being that the presence of sur¬ 
rounding matter is ignored. 

The discussion in the previous paragraph assumes many properties of so¬ 
lutions to the field equations not yet proven with mathematical rigor. Chief 
among them are that cosmic censorship holds in these scenarios, and that 
any black hole that forms in our universe (specifically here via the merger of 
two black holes, but implicitly also by processes that led to the initial black 
holes) evolves to a geometry locally describable by a unique member of the 
Kerr family (again modulo perturbations from the exterior universe). Other 
than intrinsic theoretical interest to understand merger geometries, finding 
numerical solutions for specific examples can provide strong evidence for 
these assumptions. However, the most pressing reason to study the binary 
black hole problem in recent years has been to support the effort to observe 
the universe in the gravitational wave spectrum. As discussed in more detail 
elsewhere in this volume (see Chapter 6), theoretical models of expected 
waveforms are necessary for successful detection and to decipher the prop¬ 
erties of sources. A host of tools have been developed to tackle this problem 
for black hole mergers, including post Newtonian expansions, black hole per¬ 
turbation theory, the effective one body (EOB) approach, and the geodesic 
self-force problem applicable to extreme mass ratio mergers. For comparable 
mass ratio mergers, perturbative methods break down near coalescence, and 
this is where numerical relativity contributes most to the problem. The rest 
of this section is devoted to an overview of what has been learned about 
these final stages of the merger from numerical solutions, restricting to the 
four spacetime-dimension case. For more detailed reviews see [lattM]. 

One of the results that was immediately obvious from the first full merger 
simulations of equal mass, non-spinning black holes unittiKii], and since 
then for the large swath of parameter space simulated (see for example |137) ), 
is the relative simplicity of the structure in the emitted waves during the 
transition from inspiral to ringdown (see the left panel of Fig. 7.2). This is 
the regime of evolution where the strongest-held dynamics is manifest, and 
the perturbative approaches applicable before and after should be least reli¬ 
able. Certainly the perturbative inspiral calculations do break down evolv¬ 
ing forwards to merger, and similarly for extending the quasi-normal ring- 
down backwards to this time. However, there is no signihcant intermediate 
regime of dynamics between the two, and with guidance from the numeri¬ 
cal simulations, the perturbative waveforms can be stitched together with 
relatively simple matching conditions (this is a rapidly advancing sub-held; 
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Figure 7.2 Depictions of the gravitational waves emitted during the merger 
of two equal mass (approximately) non-spinning black holes |141) . Left: 
The plus-polarized component h+ of the wave measured along the axis 
perpendicular to the orbital plane. tcAH on the horizontal axis is the time 
a common apparent horizon is first detected. Right: A color-map of the 
real component of the Newman-Penrose scalar 4/4 (proportional to the 
second time derivative of /i+ far from the BH) multiplied by r along a 
slice through the orbital plane (grey is 0, toward white (black) positive 
(negative)). From top left to bottom right the time {t — tcAH)/M of each 
panel is approximately —150, —75,0, 75. 


see |1381113911140j for a few recent examples at the time this chapter was 
written). 

With regards to the issues of theoretical interest discussed above, no sim¬ 
ulation has shown a violation of cosmic censorship, and the final state, to 
within the accuracy of the simulations and the level that researchers have 
scrutinized the geometry, is a member of the Kerr family. Moreover, though 
it is unlikely that the quasi-normal mode spectrum of Kerr is able to de¬ 
scribe all possible perturbations, in cases studied to date the post-merger 
waveforms can indeed be well approximated as a sum of quasi-normal modes. 
Of course, here we have a rather restrictive class of astrophysically minded 
“initial conditions” for the perturbed Kerr black hole formed by the merger 
of two black holes. We note that a couple of studies of single black holes per¬ 
turbed by gravitational waves have also been studied numerically beyond the 
linear regime, and similar conclusions hold |142( 1143] . 




















24 


Probing Strong Field Gravity Through Numerical Simulations 


Some of the more important numbers that have been provided by numer¬ 
ical simulations include the total energy and angular momentum radiated 
during merger (and consequently the final mass and spin of the remnant 
black hole), the spectra of quasi-normal modes excited, and the recoil, or 
“kick” velocity of the hnal black hole to balance net linear momentum ra¬ 
diated. It is beyond the scope of this chapter to list all these numbers. 
However in brief, for a baseline reference, it has been found that two equal 
mass, non-spinning black holes beginning on a zero eccentricity orbit at 
“infinite” separation radiate ~ 4.8% of the net gravitational energy dur¬ 
ing inspiral, merger and ringdown, ultimately becoming a Kerr black hole 
with dimensionless spin parameter a ~ 0.69 (due to the symmetry of this 
system, there is zero recoil). The waveform spectrum is dominated by the 
quadrupole mode in a spin-weight 2 spheroidal harmonic mode decompo¬ 
sition; the next-to-leading order is the octupole mode, which is strongly 
sub-dominant, though it briefly grows to an amplitude around l/5th that 
of the quadrupole mode near merger |141| (the energy of a mode scales as 
its amplitude squared). Changing the mass ratio decreases the energy radi¬ 
ated by roughly the square of the symmetric mass ratio r], the final black 
hole spin drops linearly with rj, new multipole moments in the waveform 
are introduced (reflecting the quadrupole moment of the effective energy 
distribution of the two particle source), and can produce recoil velocities as 
high as ~ 175km/s |14411145t I146t ITT7] . Introducing spin for the initial black 
holes can alter the radiated energies by up to a factor of roughly 2 (higher 
for spins aligned with the orbital angular momenta, lower otherwise) |148| . 
increase (decrease) the final spin for initial spin aligned (anti-aligned) with 
the orbital angular momentum (the largest aligned spin cased simulated to 
date begins with equal initial spins of a ~ 0.97, merging to a black hole 
with a ~ 0.95 |149j i. introduces precession of the orbital plane which cor¬ 
respondingly modulates the multipole structure of the waveform observed 
along a given line of sight [IbOj I15H I152j , and perhaps most remarkably can 
produce recoil velocities of several thousand km/s for appropriately aligned 
high-magnitude spins |15311154j . Figure 7.3 illustrates some of the results 
obtained for equal mass, fast spinning binary black holes. 

There are many astrophysical consequences of large recoil velocities, in 
particular for supermassive black hole mergers; we briefly mention a few 
here, together with some broader consequences of mergers on surrounding 
matter (for recent more detailed reviews see |15611157) 1. First, the veloci¬ 
ties for near equal mass, high spin mergers are large enough to significantly 
displace the remnant from the galactic core, or for the highest velocities 
even eject the black hole from the host galaxy altogether. This may be in 
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Figure 7.3 Recoil velocities from equal mass, spinning binary black hole 
merger simulations (circles) together with analytical fitting functions. Each 
black hole has the same spin magnitude a, equal but opposite components 
of the spin vector within the orbital plane, and 6 is the initial angle be¬ 
tween each spin vector and the orbital angular momentum. The dashed line 
corresponds to a fitting formula that depends linearly on the spins, while 
solid lines add non-linear spin contributions (from |155p . 


some tension with observations that seem to suggest that all sufficiently 
massive galaxies harbour supermassive black holes. If the system has a cir- 
cumbinary accretion disk, the recoil would carry the inner part of the disk 
with it, and this could be observable in Doppler-displaced emission lines 
relative to the galactic rest frame |158) . The near-impulsive perturbation 
to the gravitational potential in the outer parts of the accretion disk could 
lead to the formation of strong shocks, producing observable electromagnetic 
emission on timescales of a month to a year afterwards |159j . (Note that re¬ 
gardless of the recoil, the entire accretion disk will experience an impulsive 
change in potential due to the near instantaneous loss of energy from grav¬ 
itational wave emission at merger, also producing electromagnetic emission 
post-merger |16nj ). Earlier studies have suggested that prior to merger the 
accretion rate, and hence the luminosity of the nucleus, would be low as the 
relatively slow migration of the inner edge of the accretion disk decouples 
from the rapidly shrinking orbit of the binary. Post merger then, AGN-like 
emission could be re-ignited once the inner edge of the disk reaches the new 
innermost stable circular orbit (ISCO) of the remnant black hole. This will 
be displaced from the galactic center if a large recoil occurred, and could be 
observable in nearby galaxies (see for example m)- However, more recent 
simulations of circumbinary disks using ideal magnetohydrodynamics for the 
matter shows that complete decoupling does not occur, and relatively high 
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Figure 7.4 Left: Rest-mass density induced by a supermassive black hole 
binary interacting with a magnetized disk prior to when the binary “de¬ 
couples” from the disk, namely when the gravitational wave backreaction 
timescale becomes smaller than the viscous timescale (from [163) 1. Right: 
Poynting flux produced by the interaction of an orbiting binary black hole 
interacting with a surrounding magnetosphere. The “braided” jet structure 
is induced by the orbital motion of the black holes (from [165) 1. 


accretion rates can be maintained all the way to merger [16211163] . (The left 
panel of Fig. [774] illustrates a binary black hole system accreting surrounding 
gas). The binary orbit can cause a modulation in the induced luminosity of 
the system, which may be observable. A displaced central black hole will 
also have its loss-cone refilled, increasing the frequency of close encounters 
with stars and their subsequent tidal disruption by the black hole, with rates 
as high as 0.1/yr; the disruption could produce observable electromagnetic 
emission m . Yet another exciting prospect for an electromagnetic counter¬ 
part is an analog of the standard Blandford-Znajek mechanism (to extract 
rotational energy from a spinning black hole) induced by a tightening bi¬ 
nary within a circumbinary disk. In particular, numerical simulations have 
uncovered that binary black holes can extract both rotational and transla¬ 
tional kinetic energies when there is surrounding plasma [165) . This not only 
can power strong dual Poynting jets (emanating from each black hole), but 
the jets will increase in strength until merger, making them indirect “space- 
time tracers”—the right panel of Fig. 7.4 depicts the resulting “braided” 
structure of the Poynting flux. 

As a final comment we note that the majority of work, both numerical and 
analytic, has been devoted to studying zero-eccentricity mergers, due to the 
prevailing view that these will dominate event rates. However, there are bi¬ 
nary formation mechanisms that can produce high-eccentricity mergers (see 
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the discussion in |166j for an overview and references). One of the interest¬ 
ing results from the handful of studies including large eccentricity performed 
to date HSHIMUM] is that zoom-whirl orbital dynamics is possible for 
comparable mass binaries. In the test particle limit, zoom-whirl orbits are 
perturbations of the class of unstable circular geodesics that exist within the 
ISCO; further, they exhibit extreme sensitivity to initial conditions where 
sufficiently fine-tuned data can exhibit an arbitrary number of near-circular 
“whirls” at periapse for a fixed eccentricity geodesic. Away from the test 
particle limit gravitational wave emission adds dissipation to the system, 
though what the simulations show is that even in the comparable mass limit 
the dissipation is not strong enough to eradicate zoom-whirl dynamics, but 
merely limits how long it can persist. 


7.3.3 Black Hole-Neutron Star/Binary Neutron Star Mergers 

Non-vacuum compact binary systems-i.e., those involving at least one neu¬ 
tron star-are also the subject of intense scrutiny. These systems produce 
powerful gravitational waves and likely also lead to intense neutrino and 
electromagnetic emission that could be detected by transient surveys or by 
dedicated follow up by the astronomical community. In particular they are 
posited to be the progenitors of short gamma ray bursts (sGRBs) and a 
host of other transient phenomena [Hoi in]. Signals from these systems 
can thus carry a wealth of information about gravity, the behavior of mat¬ 
ter at nuclear densities, and binary populations and their environments. The 
challenge for simulations is to obtain predictions to confront with observa¬ 
tions. 

Relative to the two black hole case, the most obvious complication in 
the simulation of binaries with neutron stars is the need to include non- 
gravitational physics. The simplest relativistic model of a neutron star cou¬ 
ples relativistic hydrodynamics to the Einstein equations and, using a sim¬ 
plified equation of state (EOS), the first successful simulations of binary 
neutron star mergers within this framework were presented in [1721 1173) . 
Since the time of those studies, the community has made steady progress in 
exploring the full parameter space relevant to astrophysical mergers, while 
simultaneously increasing the fidelity of the matter modeling through inclu¬ 
sion of the electromagnetic interaction, neutrino and radiation transport, 
nuclear reactions, and other physics. A crucial unknown here is the EOS 
that describes matter at nuclear densities: it plays a leading role in the phe¬ 
nomenology of the system as, for a given stellar mass, it regulates the star’s 
radius, affects its response to tidal forces, and affects its ability to resist 
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collapse to a black hole when it accretes matter (or collides with another 
star). Given the difficulty of first-principles calculations or probing similar 
conditions in laboratories, detailed knowledge of the nuclear density EOS 
is likely to come only through astronomical observations, and prospects for 
doing this through gravitational waves are particularly exciting—see for ex¬ 
ample [17411175) . 

While the pericenter is large these systems evolve much like black hole bi¬ 
naries. The orbit shrinks due to the emission of gravitational radiation with 
internal details of the stars playing essentially no role. However, finite body 
effects become important as the orbit tightens. In the remainder of this para¬ 
graph we focus on binary neutron stars, returning to black hole-neutron star 
systems in the following paragraph. Tidal forces deform both stars (which 
can even induce crust-shattering |176| ). leaving subtle imprints in the ensu¬ 
ing gravitational waves. This behavior intensifies until the point of merger, 
when the local velocities reach a sizable fraction of the speed of light, ending 
in a violent collision that ejects neutron rich matter due to shock heating and 
extreme tidal forces. Figure [7^ (left panel) illustrates waveforms obtained in 
an equal mass binary neutron star system for different EOS models, demon¬ 
strating how significantly this can affect the behavior. In general terms, for 
the typically expected neutron star masses of 1.2 — 1.8Mq the merger yields 
a hot, differentially rotating, hyper-massive neutron star (HMNS). Such an 
object will promptly collapse to a black hole if the total binary mass is above 
2.6 — 2.8Mq, depending on the stiffness of the EOS. Otherwise, a delayed 
collapse takes places as the star is initially supported by differential rota¬ 
tion and thermal pressure. During this stage, the HMNS rotates and emits 
gravitational waves with frequencies in the range 2 < / < 4Khz, with a 
characteristic frequency proportional (and relatively close) to the Keplerian 
velocity (MhmnsZ-^hmns)^^^ (s-S-j [IZIj)- ^ scale of tens of milliseconds 
however, such support diminishes due to gravitational radiation, angular 
momentum transport via hydrodynamical and electromagnetic effects, and 
cooling due to neutrino emission (these effects have just begun to be studied, 
e.g. |1781117911180) 1. The black holes that form in both prompt or delayed 
cases are (reasonably) well-described by a Kerr solution with a spin pa¬ 
rameter J/M^ < 0.8, surrounded by left-over material, much of which is 
bound (e.g., [181) 1 and can form an accretion disk with a mass on the order 
of ~ 0 — 0.3Mq. The amount depends on the EOS, mass ratio, and elec¬ 
tromagnetic fields (though this latter effect is still largely unexplored) and 
is intuitively larger for longer lived HMNS as more angular momentum is 
transferred outwards to the envelope. Importantly, this is enough material 
to form a sufficiently massive disk as called for in models of sGRBs. Some 
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material will be ejected (again, the amount depending upon various param¬ 
eters) and will decompress to form heavy elements through the r-process; if 
these merger events are frequent this could account for a significant fraction, 
if not the majority of such elements in the Universe. Subsequent decay of 
the more radioactive isotopes could lead to a so-called kilo- or macronova 
(reports of the afterglow of the recent sGRB 130603B are consistent with 
this |182l I183| ). Observation of these signatures together with gravitational 
wave observations will allow us to make contact between simulations and 
the birth of a black hole. However, gravitational waves emitted during the 
HMNS and collapse stages have a higher frequency than those reachable by 
LIGO/VIRGO/KAGRA, and will take third-generation facilities to detect. 
Nevertheless, up to the frequencies that existing (and near future) detectors 
can probe, subtle differences in the gravitational waveforms should allow 
for constraining the radius of the neutron stars to within 10% |184| 1185) . 
Simulations are further probing possible counterpart signals from neutrino 
production |186j and electromagnetic interactions [1871118811189] . 

Black hole-neutron star binaries display even more complex merger dy¬ 
namics. Indeed, at an intuitive level one expects significant differences to 
arise depending on whether the tidal radius Rt (oc Rns (SMbh/ATns)^^^) 
lies inside or outside the black hole’s inner most stable circular orbit ra¬ 
dius (Risco)) which ranges from Mbh to 9 Mbh for a prograde to retrograde 
orbit about a maximally spinning black hole. This is clearly borne out in 
simulations exploring a range of mass ratios and black hole spins, show¬ 
ing markedly different behavior in the ensuing dynamics and gravitational 
waves produced. Qualitatively, for sufficiently high spins and/or sufficiently 
low mass ratios, the star signihcantly disrupts instead of plunging into the 
black hole. As a result, gravitational waves promptly “shut-off” at a fre¬ 
quency related to the star’s EOS. Figure 7.5 (right panel) illustrates this for 
different EOS models in a 3 : 1 mass ratio black hole-neutron star system. 
When disruption occurs during the merger, a significant amount of mate¬ 
rial, in the range 0.01 — 0.3Mq, can remain outside Risco- This material 
will be on trajectories having a range of eccentricities, with the fraction that 
is bound falling back to accrete onto the black hole at a rate governed by 
the familiar law M oc |192t 1193] . The details however depend on many 
factors, including spin-orbit precession as illustrated in Fig. 7.6 The matter 
that is ejected (< O.OSMq) can be have speeds up to ~ 0.2c [19411191) . This, 
together with the amount of likely accretion, is in the range assumed by 
models predicting that black hole-neutron stars mergers can power sGRBs, 
kilonovae and related electromagnetic counterparts. Consequently, a simi¬ 
lar array of electromagnetic signatures and r-process elements could result 
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Figure 7.5 Examples of the “plus” polarization component of gravitational 
waves from binary neutron star mergers, measured 100 Mpc from the 
source along the direction of the orbital angular momentum. The differ¬ 
ent curves correspond to different choices of the EOS of the neutron star 
matter, labeled APR4, ALF2, H4 and MSI. For a IAMq neutron star, 
the APR4, ALF2, H4, MSI EOS give radii of 11.1,12.4,13.6,14.4km re¬ 
spectively. Left: Mergers of an equal mass binary neutron star system 
(with mi = m 2 = I. 4 M 0 ). A hypermassive neutron star (HMNS) is 
formed at merger, but how long it survives before collapse to a black 
hole strongly depends on the EOS. The H4 case collapses to a black hole 
« 10ms after merger; the APR and MSI cases have not yet collapsed 
~ 35ms after merger when the simulations where stopped (the MSI EOS 
allows a maximum total mass of 2 . 8 M 0 , so this remnant may be sta¬ 
ble). The striking difference in gravitational wave signatures is self evident 
(from |190j L Right: Emission from black hole-neutron star mergers, with 
iTT-Bii = 4.O5M0,mNs = 1.35M0. Variation with EOS is primarily due to 
coalescence taking place earlier for larger radii neutron stars (from [191j l. 


as with binary neutron star mergers, and the gravitational wave signals 
could be ideal to differentiate between them. For the subset of black hole 
neutron star mergers where Rt < .RiscO; the star plunges into the black 
hole with little or no material left behind, and the resulting gravitational 
wave signal will be much like that of a binary black hole system with the 
same binary parameters. Counterparts such as sGRBs or kilonova requiring 
significant accretion disks or unbound matter are therefore not favored for 
this sub-class of binary. Nevertheless, interesting electromagnetic precursors 
could be induced by magnetosphere-black hole interactions prior to merger 
(e.g. 1193 [Mira). 

An important observation is that black hole-neutron star systems are, in 
all likelihood, more massive than binary neutron star systems. Therefore, the 
wave frequency peaks at lower frequencies than binary neutron star merg¬ 
ers, offering better prospects for observing non-linear effects by near-future 
detectors. Indeed, since the characteristics of gravitational waves depend on 
masses, spins and the EOS, black hole-neutron star systems provide perhaps 
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Figure 7.6 General relativistic hydrodynamic simulations of the merger of 
a 9.8Mq, a = 0.9 black hole with a 1.4M0 neutron star, from |193) . The top 
two panels are from a case where the spin and orbital angular momentum 
vectors are aligned; the bottom two where the initial (~ 9 orbits before 
merger) misalignment is 40°. The left two panels are at a time when half 
the material has been absorbed by the black hole, showing matter densities 
above ^ 6 x 10^°g/cm^; the right two are 5ms later, showing densities 
above ^ 6 x 10®g/cm^, and the facing quadrant has been cut from the top- 
right rendering. These results illustrate the profound affect spin-induced 
precession can have on the matter disruption and subsequent accretion. 


the best prospects for extracting key physical information about neutron 
stars |175t 1198] . To date the majority of simulations have focused on a black 
hole spin aligned with the orbital angular momentum, with the exception 
of |193j . which showed that the above conclusions hold qualitatively even 
with inclinations < 30° of the spin axis away from alignment. For larger 
inclinations, of the disrupted material a smaller fraction forms a disk on a 
short timescale following merger, while a larger fraction follows an eccentric 
trajectory and returns to interact with the black hole on longer timescales. 

As in the binary black hole case, incipient efforts are examining encounters 
with high initial orbital eccentricity in non-vacuum binaries (e.g. [19911200] ). 
Qualitatively, much of the same phenomenology of outcomes can occur as 
with quasi-circular inspirals (except that now zoom-whirl orbital dynamics is 
also possible), though the details can be drastically different. For example, 
in high eccentricity encounters of a neutron star with a black hole, tidal 
disruption can occur for higher mass ratio systems and smaller black hole 
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spin, as the effective inner most stable orbit is closer to the black hole for 
eccentric orbits. There can also be multiple, partial disruptions on each of 
the last several periapse passages, ejecting larger amounts of material and 
leaving behind more massive accretion disks than otherwise possible. On 
close periapse passages (even without disruption) f-modes can be impulsively 
excited in the star, or both stars in a binary neutron star encounter. These 
modes are too low amplitude/high-frequency to be directly observed with 
the current generation of ground-based gravitational wave detectors, though 
they may indirectly be measured in the leading order part of the waveform, 
since from the perspective of the binary the f-modes are a new channel of 
energy dissipation. The impulsive tidal interaction may also cause crust- 
shattering m, leading to electromagnetic emission similar to the resonant 
excitation induced shattering in quasi-circular inspiral |176j . 

For merger simulations involving neutron stars, the current frontier is to 
add more matter physics to the models (resistive magnetohydrodynamics, 
radiation and neutrino physics, multi-component fluids, “realistic” high tem¬ 
perature equations of state, etc.). Given the many orders of magnitude of 
spatial and temporal scales involved, as well as the complexity of the micro¬ 
physics, it will likely be several years before both realistic models and the 
computational power necessary to simulate them accurately are available. 
Due to space constraints we will not list all the directions currently being 
pursued, referring the reader to recent reviews in |2n2( 1203] 1204] . 


7 . 3.4 Gravitational Collapse to a Neutron Star or Black Hole 

Considerable efforts have been undertaken to study gravitational collapse 
to a neutron star or a black hole, in particular within the context of core¬ 
collapse supernovae. Here, stars with masses in the range IOMq < M < 
IOOMq at zero-age main sequence form cores which can exceed the Chan¬ 
drasekhar mass and become gravitationally unstable. This leads to collapse 
which compresses the inner core to nuclear densities, at which point the full 
consequences of general relativity must be accounted for. Depending upon 


the mass of the core, it can “bounce” or collapse to a black hole. Figure 7.7 


displays representative snapshots of the behavior of a collapsing 75Mq star 
at different times. The collapse forms a proto-neutron star which later col¬ 
lapses to a black hole. In the case of a bounce, an outward propagating 
shock wave is launched which collides with still infalling material and stalls. 
Observations of core-collapse supernovae imply that some mechanism is ca¬ 
pable of reviving the shock, which is then able to plow through the stellar 
envelope and blow up the star. This process is extremely energetic, releasing 
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energies on the order of lO^^erg, the majority of which is emitted in neu¬ 
trinos. For several decades now, the primary motivation driving theoretical 
and numerical studies has been to understand what process (or combination 
of processes) mediates such revival, and how (for a recent review see [205) ). 
Several suspects have been identified; heating by neutrinos, (multidimen¬ 
sional) hydrodynamical instabilities, magnetic fields and nuclear burning 
(see e.g. [20611207] i. With the very disparate time and space scales involved, 
a multitude of physically relevant effects to consider, and the intrinsic cost to 
accurately model them (e.g., radiation transport is a 7-dimensional problem) 
progress has been slow. Moreover, electromagnetic observations do not pro¬ 
vide much guidance to constrain possible mechanisms, as they can not peer 
deep into the central engine. On the other hand, observations of gravitational 
waves and neutrinos have the potential to do so, provided the explosion is 
sufficiently close to us. Thus, in addition to exploring mechanisms capable of 
reviving the stalled shock, simulations have also concentrated on predicting 
specific gravitational wave and neutrino signatures. 


Modeling gravity using full general relativity has only recently been un¬ 
dertaken |2U8| . though prior to this some of the more relevant relativistic 
effects were incorporated (e.g. [MUSiniEIIlEIS])- While the full resolution 
of the problem is still likely years ahead, interesting insights into fundamen¬ 
tal questions and observational prospects have been garnered. For example, 
simulations have shown that in rotating core collapse scenarios, gravitational 
waves can be produced and their characteristics are strongly dependent on 
properties of the collapse: the precollapse central angular velocity, the de¬ 
velopment of non-axisymmetric rotational instabilities, postbounce convec¬ 
tive overturn, the standing accretion shock instability (SASI), protoneutron 
star pulsations, etc. If a black hole forms, gravitational wave emission is 
mainly determined by the quasi-normal modes of the newly formed black 
hole. The typical frequencies of gravitational radiation can lie in the range 
~ 100 — 1500Hz, and so are potential sources for advanced earth-based grav¬ 
itational wave detectors (though the amplitudes are sufficiently small that it 
would need to be a galactic event). As mentioned, the characteristics of these 
waveforms depend on the details of the collapse, and hence could allow us to 
distinguish the mechanism inducing the explosion. Neutrino signals have also 
been calculated, revealing possible correlations between oscillations of grav¬ 
itational waves and variations in neutrino luminosities. However, current 
estimates suggest neutrino detections would be difficult for events taking 
place farther than kpc distances [208] . 
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Figure 7.7 Density colormaps of the meridional plane of a collapsing ToMq 
star superposed with velocity vectors at various times after bounce (and 
with different spatial ranges to zoom-in on particularly relevant behavior). 
The collapse first forms a proto-neutron star which later collapses to a BH 
(shown in the bottom panels). (From |213) b 


7.3.5 Ultra-relativistic Collisions 

Some of the early interest in the ultra-relativistic collision problem stemmed 
from investigations by Penrose |214j into its relevance to questions of cosmic 
censorship. It was known that collisions of gravitational waves with planar 
symmetry in 4d lead to the formation of naked singularities regardless of 
how “weak” the initial curvature. This is not considered a serious counter¬ 
example to cosmic censorship as the spacetime is not asymptotically flat, nor 
for that matter are there black hole solutions with planar symmetry in 4d 
vacuum Einstein gravity (with zero cosmological constant), so in a sense the 
question of censorship is not particularly meaningful here. However, taking 
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the infinite boost limit of the Schwarzschild metric (scaling the rest mass 
m to zero as the boost 7 —?• oo while the energy E = m'j remains finite 
ms]) results in the Lorentz contraction of the curvature to a plane-fronted 
gravitational shock wave, with Minkowski spacetime on either side. One can 
then consider what happens when two such shock waves, traveling in op¬ 
posite directions, collide. Given the resemblance between the two scenarios, 
the infinitely boosted black hole collision is a natural place to test cosmic 
censorship, especially since the geometry approaches Minkowski spacetime 
transverse to the center of each shock sufficiently rapidly to remove the 
trivial objections to the plane-symmetric gravitational wave collisions. Pen¬ 
rose found a trapped surface in a zero-impact parameter, infinite 7 black 
hole collision, and though the metric to the causal future of the collision is 
unknown, this is a good indication that cosmic censorship does holds here. 

More recently two additional lines of research have come to the fore moti¬ 
vating the study of ultra-relativistic collision geometries. The first is a con¬ 
sequence of the observation that if extra spatial dimensions exist, then the 
true Planck scale could be much different from the effective 4-dimensional 
scale one would otherwise expect [21611217] . In particular, a “natural” solu¬ 
tion to the hierarchy problem results if the Planck energy is on the order of 
a TeV. If that is the case it was conjectured that particle collisions at the 
Large Hadron Collider (LHC) and in cosmic ray collisions with the earth 
with center of mass energies above this could result in black hole forma¬ 
tion |218t I219j The conjecture is essentially based on two premises: that 
Thorne’s hoop conjecture |223j can be applied to the collision to deduce 
whether the purely classical gravitational interaction between particles will 
cause a black hole to form, and if so, that the quantum interactions are 
sufficiently “local” to not alter this conclusion (until Hawking evaporation 
becomes significant). The second motivation comes from applications of the 
AdS/CFT correspondence of string theory to attempt to explain the for¬ 
mation and early time dynamics (before hadronization) of the quark-gluon 
plasma formed in relativistic heavy ion collisions (RHIC) (see Sec. 7.3.6 for 
more on this). Here, the gravitational dual to a heavy ion collision is conjec¬ 
tured to be an ultrarelativistic black hole collision in the bulk asymptotically 
AdS spacetime. 

The first 4d ultra-relativistic head-on black hole collision simulations (up 
to 7 ~ 3) were carried out in |224j . followed by several studies with general 


To date, the LHC has not seen evidence for black hole formation in searches of collisions with 
center of mass energies up to 8 TeV |22QI 1221) : likewise, no signs of black hole formation have 
yet been observed in cosmic ray collisions 1^ . the most energetic of which can have much 
larger center of mass energies. 
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impact parameters [22511226j , including the effects of black hole spin [2271 
1228] , and collisions in higher dimensions [229] . A wealth of interesting results 
have emerged, a select few of which we briefly summarize here. Outcomes 
of most interest to LHC searches include the critical impact parameters 
for black hole formation, and the energy and angular momentum lost to 
gravitational waves as a function of impact parameter. This determines the 
formation cross-section and initial spectrum of black hole masses that will 
subsequently Hawking evaporate. Extrapolated results from 4d head-on col¬ 
lisions give 14 ± 3% energy emitted in gravitational waves, roughly 1/2 the 
Penrose trapped surface calculation, though consistent with the 16% ob¬ 
tained using perturbative analytic methods [230]. As the impact parameter 
increases, the radiated energy and now angular momentum increases, though 
the former is still less than trapped surface estimates [231| . Qualitative fea¬ 
tures of the spectrum of emitted waves can be understood appealing to 
the analytic zero-frequency and point particle limit calculations |232| . The 
largest gravitational wave fluxes arise near the threshold impact parame¬ 
ter. Here, in the 4d case, the binary exhibits behavior akin to zoom-whirl 
dynamics of black hole geodesics, though not in 5d (presumably due to the 
stronger effective gravitational potential, related to the fact that there are no 
stable circular orbits about Myers-Perry black holes in dimensions greater 
than 4) [22611229] . 

Because of the zoom-whirl behavior in 4d, it was argued in [167] that turn¬ 
ing to the threshold in the large 7 limit essentially all the initial kinetic en¬ 
ergy of the black holes would be converted and radiated out as gravitational 
waves. However, the results presented in |228j show this is likely not true, 
due to what appears to be strong self-absorption of the emitted gravitational 
energy by the black holes. These simulations only went to 7 ~ 2.5; however 
if they in fact provide a decent approximation of the large 7 limit,then one 
concludes that as much as 1/2 the kinetic energy could be converted to rest- 
mass energy in the black holes (the rest to gravitational waves), even in close 
scattering encounters (which is in fact consistent with a perturbative calcu¬ 
lation in the extreme mass ratio limit [233j 1. The surprising consequence of 
this is that two “microscopic” black holes each of rest-mass m scattered off 
one another with 7 ;§> 1 and finely tuned impact parameter could grow to 
two “macroscopic” black holes moving apart sub-relativistically, each with 
rest-mass ~ 7717 / 2 . 

A further intriguing result for the 5d case presented in |229j is that for a 
small range of impact parameters near threshold, curvature invariants grow 
rapidly at the center of mass shortly after what appears to be a scattering 
event, though the code crashed before the final outcome could be determined. 
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No encompassing apparent horizon is detected then, which could simply be 
because of the nature of the time coordinate employed, or could be due 
to a naked singularity that is forming. If the former, this would be a new 
outcome to the black hole scattering problem in 5d, namely three black holes; 
if the latter, this would be another example (in additional to the Gregory- 
Laflamme instability of black strings, and possible prolate dust collapse [22]) 
showing a violation of cosmic censorship. 

The high-speed limit has also shed some light on the mechanism respon¬ 
sible for large recoil velocities seen in merger simulations of astrophysically 
relevant inspirals with certain spin configurations (see the discussion in 
Sec. 7.3.2). In particular, there have been suggestions that the large re¬ 
coils require the formation of a common horizon to effect the gravitational 
wave emission of “field momentum” associated with what is otherwise purely 
kinematical properties of the orbit; in reaction the merger remnant receives 
a kick in a direction that conserves linear momentum |234j . However, in 
high speed merger simulations with similar black hole spin setups, even in 
scattering cases where a common horizon does not form, large recoils are 
observed |227j . This is consistent with the heuristic explanation of the su¬ 
perkicks presented in m as arising from frame-dragging induced Doppler 
boosting of the radiation emitted by the binary motion. 

The first ultra-relativistic collision simulations of “solitons” (non-singular 
compact distributions of matter) were carried out in |235j , consisting of the 
head-on collision of two boson stars each with compactness 2M/R ~ 1/20, 
and center of mass boosts up to 7 = 4 (u ~ 0.968). The main goal of 
the study was to test the hoop conjecture arguments for black hole for¬ 
mation; hence the use of boson stars as model particles given that their 
self-interaction is weak compared to gravity in this limit. Black hole forma¬ 
tion was observed above a critical boost ss 3, roughly one third the value 
predicted by the hoop conjecture. Similar results were later obtained using 
an ideal fluid (fermion) star as the model particle [23611237| . The study in 
|236| used less compact stars that pushed the critical boost to 7 c ~ 8.5, but 
again found this to be a similar factor less than the hoop conjecture esti¬ 
mate. It was argued that the lower thresholds are due to the compression 
of one particle by gravitational focusing of the near-shock geometry of the 
other particle, and vice-versa. This conclusion was anticipated by a geodesic 
model of black hole formation presented in |238j . It is remarkable that such 
a simple model, and for that matter the trapped surface calculations as well, 
predict the qualitative properties of what is ostensibly the regime where the 
most dynamical, non-linear aspects of the Einstein equations are manifest. 
On the other hand, a recent calculation of the gravitational self-force us- 
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ing effective field theory techniques in the large 7 limit show that many 
simplifications arise here; in particular the non-linear interactions coming 
from gravitational bulk vertices are suppressed by factors of 1 / 7 “^ [239| (see 
also |24n| ). Aside from giving strong evidence that the hoop conjecture is 
applicable to the classical collision problem, these studies also support the 
expectation that the outcome of sufficiently supercritical 7 > 7c collisions 
are insensitive to the details of the particle self-inter actions. This is essential 
for black hole formation in super-Planck particle collisions to be a robust 
conclusion, despite the lack of detailed calculations in a full quantum (grav¬ 
ity) theory. This also justifies the use of black holes as model particles, 
which from a classical gravity perspective is (in theory) a simpler problem 
to simulate, due to the absence of matter. 

The motivation and applications of the AdS/CFT correspondence in string 


theory are discussed below in 7.3.6 here we briefly comment on what RHIC- 


motivated studies have taught us about ultra-relativistic collisions. The rel¬ 
evant spacetime for this problem is 5d AdS, and in particular the Poincare 
wedge, as its boundary is conformal to 4d Minkowski spacetime. Solving 
the full Einstein equations in 5 spacetime dimensions without symmetries 
and resolving the geometry dual to highly boosted concentrations of energy 
would be an extremely challenging computation to perform. To date then, 
existing studies (see m for a review) have made simplifying approxima¬ 
tions: each particle is modeled as a finite-width gravitational wave with pla¬ 
nar symmetry transverse to the collision axi^ This effectively reduces the 
numerical evolution to 2 -|- 1 dimensions, and characteristic approaches have 
proven highly successful for this problem. Though the topology and asymp¬ 
totics are quite different from the 4d asymptotically flat case, there is some 
similarity. Most relevant to this discussion is that the infinite boost limit is 
similar to the Penrose/Aichelburg-Sexl superposed shock-wave construction; 
in both cases trapped surfaces can be found |242j . yet the full solution to 
the causal future of the collision is unknown. The numerics have solved the 
finite width planar collision problem, showing that a black hole (with planar 
topology) does form in this case, and resolving the spacetime to the future of 
the shock. In particular, post-collision along the future lightcone of the colli¬ 
sion, the amplitude of the shock, as projected onto the Minkowski boundary, 
decays as a power law in time; within the lightcone, after a time roughly 
consistent with inferred thermalization times in RHIC experiments, the near 
boundary metric fluctuations transform to a state that can be characterized 


® Note added in proof : a first simulation without symmetry assumptions in 5d AdS was 
presented in Em 
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as an expanding, cooling hydrodynamic flow [76ll243j . For further details on 
numerical relativity applications in the realm of high energy see [ 213 . 


7 . 3.6 Gravity in d ^ 4, 

Beyond ultrarelativistic collisions, numerical relativity has also been crucial 
in exploring the behavior of gravity in both stationary and time dependent 
scenarios beyond d = A. There are several motivations for doing so. On 
one end there is the desire to understand gravity at a fundamental level by 
contrasting known behavior in d = 4 to what arises in different dimensions. 
Higher dimensions are required by string theory, and this has inspired many 
speculative theories: for example TeV scale gravity/braneworIds [24511217] . 
some models of inflation |246j and modern cyclic models |247) of the universe. 
Lower dimensions have also been used to provide a simpler setting to gain 
intuition about quantum gravity (e.g., in 2 + 1 [248j and in 1 + 1 dimensional 
dilaton gravity |249| 1. At the other end, compelling practical reasons are 
provided by the role gravity may play to understand phenomena described 
by field theories through holography |250t I251j . 

For more information, readers are directed to the recent book |252j . Here, 
for brevity we mainly focus on time-dependent problems, though we briefly 
review stationary solutions; in particular, those that are relevant to existing 
or future dynamical studies. 

Black holes in dimensions d> A 

Understanding the landscape of stationary solutions with event horizons has 
been the focus of considerable effort [2531 1254] . This work has illustrated 
how much richer the space of stationary black object solutions in higher 
dimensions is compared to the d = A case. A case in point is the broader class 
of topological structures allowed, which includes hyper-spherical black holes, 
black rings and a combination of these latter two giving “black Saturns”, 
black strings, black branes, etc. Interestingly, several topologically distinct 
solutions can have the same asymptotic charges, showing some degree of 
non-uniqueness of black hole solutions in higher dimensions. However, a 
particularly intriguing conjecture is that uniqueness can be restored by the 
additional requirement of stability. This possibility is implied by the fact that 
linear perturbations of many of these solutions are unstable. As a further 
contrast with d = A stationary black holes, there is no “Kerr-like” bound 
for spinning black holes in d > 5 as they can have arbitrarily large angular 
momenta. Again, related to the uniqueness issue, these ultra-spinning black 
holes are unstable [255] • Numerical solutions are required to understand the 
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Figure 7.8 A sequence of snapshots showing the evolution (left to right) of 
an unstable black string; see |256j for details. 


non-linear dynamics of unstable black objects, and to-date this has only 
been achieved for black strings in d = 5 |256j and ultra-spinning black holes 
in d G 5..8 [^1^ . 

Black strings are black hole solutions extended along a trivial (option¬ 
ally) compactified extra-dimension. For simplicity, and because it is the one 
studied numerically, we restrict to d = 5, and so the static black string is 
given by the d = 4 Schwarzschild solution cross a circle with (asymptotic) 
length L. Gregory and Laflamme showed that linearized perturbations of 
such a black string admits exponentially growing modes above some critical 
L/M (with M the mass per unit length of the black string) |259j . Further, 
thermodynamical arguments suggested that above this ratio the entropically 
preferred solution would be a d = 5 Schwarzschild-Tangherlini black hole. 
Thus, it appeared possible that the effect of these growing modes would be 
to eventually cause the black string to pinch-off and give rise to a spheri¬ 
cal configuration. Naturally, if that happened, cosmic censorship would be 
violated, indicating yet again that gravity in d = 4 is rather special. 

To understand the dynamical behavior of the solution, a full non-linear— 
and so necessarily numerical—analysis is required. Such a study was pre¬ 
sented in [256], and revealed that the instability unfolds in a self-similar 
fashion, where the black string horizon at any given time could be seen as 
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thin strings connected by hyperspherical black holes of different radii (see 
Fig [7^. As the evolution proceeds, pieces of the string shrink further while 
others give rise to spherical black holes bulges, and the horizon develops a 
fractal structure. Interestingly, such behavior is reminiscent of the one dis¬ 
played by a thin column of fluid through the Rayleigh-Plateau instability (see 
[260| '). In the case of the black string, extrapolating the numerical results 
shows that the ever-thinning string regions eventually reach zero size, reveal¬ 
ing a massless naked singularity in finite time. Thus, perturbed black strings 
do provide a counter-example to the cosmic censorship conjecture, though in 
d = 5. In still higher dimensions, the outcome is expected to be qualitatively 
similar up to a critical dimension beyond which stable, non-uniform black 
strings states are entropically favored. Perturbative analysis indicates that 
the critical dimension is d = 13 ISU, though recent work making use of a 
local Penrose inequality suggests it may be as low as d = 11 |262j . 

This result has application beyond black string spacetimes, as many of 
the higher dimensional black hole solutions have a near horizon geometry 
that can be mapped, in appropriate regions, to black strings. For instance, 
ultra-spinning black holes satisfy the Gregory-Laflamme instability condi¬ 
tion around the polar region |255j . Such black holes are thus expected to 
develop growing deformations about the poles of the horizon when per¬ 
turbed. These can be both axisymmetric modes that would evolve toward ax¬ 
ially “pinched” or ring-like configurations |263) . and also non-axisymmetric 
modes. The latter however would induce a time-varying quadrupole moment 
that would radiate angular momentum, allowing for the possibility that grav¬ 
itational wave emission could regulate the instability, in particular if the 
non-axisymmetric modes are the dominant unstable ones. Numerical simu¬ 
lations for systems with angular momentum mildly higher than the critical 
value show this is precisely the case [Ml ESS], where a “bar-mode” devel¬ 
ops that radiates angular momentum until the black hole settles down to a 
sub-critical, stable state. For larger initial spins it has been speculated that 
non-axisymmetric modes can grow more rapidly than gravitational wave 
emission can reduce the spin to sub-critical, and the horizon might then 
fragment into multiple pieces [25511258j . These cases have yet to be explored 
beyond the linear level. 

As a last example we mention that solutions describing large black holes 
in Randall-Sundrum models were numerically constructed in |264j . disprov¬ 
ing a conjecture that such solutions could not exist [2651 I266j . Moreover, 
the particular Ricci flow method employed to obtain the solutions argues 
implicitly in favour of their stability. 
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AdS/CFT duality applications 


The AdS/CFT correspondence |25nt I251| provides a remarkable framework 
to study certain strongly coupled gauge theories in d dimensions by mapping 
to weakly coupled gravitational systems in d + 1 dimensions. A large body 
of work has been built since the introduction of this correspondence; we 
will not review it here. Rather, we concentrate on a handful of applications 
where numerical simulations have been crucial to the understanding of the 
gravitational aspects of the problems. The relevant spacetimes typically in¬ 
volve black holes, and are asymptotically AdS, the latter property creating 
delicate issues on both analytical and numerical fronts. This is in part due 
to the timelike nature of the AdS boundary, with the consequence that the 
correct specification of boundary data (in addition to the initial configura¬ 
tion) is crucial for a well-defined evolution that can be mapped to the CFT 
description of the problem. The boundary conditions can be derived from 
the limiting behavior of the Einstein equations approaching the boundary, 
together with constraints from imposing that the spacetime approaches AdS 
at the appropriate rate in a suitable gauge. 


Many interesting applications have been pursued using holography, spurred 
by work beginning soon after the formulation of the correspondence in¬ 
dicating the rather spectacular breadth of possible applications to finite- 
temperature field theory (see |267) for a review). Some highlights include 
that the hydrodynamic behavior of field theory is captured by correlation 
functions in the low-momentum limit, that hydrodynamics modes in rele¬ 
vant field theory states correspond to low-lying quasinormal modes of an 
AdS black brane solution, and for such states there is a universal viscos¬ 
ity to entropy density ratio pis = I/Att for a broad class of theories with 
gravitational duals. The value of p/s is remarkably close to that inferred 
from hydrodynamic models of the quark-gluon plasma (QGP) formed in 
relativistic heavy ion collisions, and this observation has led to the new 
approach of using AdS/CFT to try to understand the QGP (for a review 
see [268) 1. Though the N = A SYM theory of the duality is not deconfined 
QCD, there are sufficient similarities that one might hope the former can 
give insights into aspects of the problem not easily calculable via traditional 
techniques (perturbative Feynman diagrams and lattice QCD). For exam¬ 
ple, using AdS/CFT purely gravitational studies can be used to estimate 
the thermalization timescales post-collision, and the subsequent evolution 
of the expanding plasma to the point of hadronization. As mentioned above 
in Sec. 7.3.5, a series of groundbreaking works |269) 1270) 1271) studied the 


behavior of the spacetime when two gravitational shock waves collide. This 
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is expected to offer a decent approximation to the dynamics along the beam 
axis in central (head-on) collisions. The results for the thermalization time 
(one definition of which is the time after collision when the boundary stress 
tensor of the CFT is well approximated by the leading order terms in the hy¬ 
drodynamic expansion) are broadly consistent with the times inferred from 
experiment. Moreover, the subsequent hydrodynamic flow exhibits a form 
of boost-invariance similar to the predicted Bjorken flow [2721127,*^ . the dif¬ 
ference being characterizable as a modest dependence of the energy scale 
of the flow on rapidity [76l I243j . Figure 7.9 illustrates the energy density 
measured at the AdS boundary from simulations on the gravitational side; 
the left image is from a shock collision simulation, and the right is the re¬ 
laxation of a highly perturbed black hole. Soon after the collision (left) and 
from the beginning of ringdown (right), a hydrodynamical description on 
the field theory side matches the observed near-boundary metric behavior 
to an excellent degree. 

Another front where the duality is being exploited is to understand the 
behavior of a system in the ground state of a given Hamiltonian when a 
“quenched interaction” is introduced. Here the response of an initial thermal 
equilibrium state of the theory under rapid variations of suitable operators 
can be studied using the correspondence. As the behavior on the gravita¬ 
tional side is governed by the dynamics of an appropriately perturbed black 
hole, a universal response is uncovered which, on the CFT side, means the 
system responds in a way only dependent on the conformal dimension of 
the quench operator in the vicinity of the ultraviolet hxed point of the the¬ 
ory ^m\Tm\Tm\ . 

As a last example we mention an application of the duality in the op¬ 
posite direction; using knowledge of the behavior on the field theory side 
to discover and analyze novel features on the gravitational side. It is well 
known that field theories at sufficiently high energies admit a hydrodynam¬ 
ical description; this motivated works that established a duality between 
gravity and hydrodynamics for relativistic, conformal fluids. Specifically, it 
was shown that in the limit of long wavelength perturbations of black holes 
the Einstein equations projected onto the AdS boundary reduce to the fa¬ 
miliar relativistic hydrodynamics for a viscous fluid (e.g. [277] 1. Numerical 
work has demonstrated that the hydrodynamic description matches the be¬ 
havior of full, non-linear solutions of the Einstein equations surprisingly well 
in may situations as mentioned above with RHIC applications (and see Fig¬ 
ure 7.9). Taken the other direction then, this duality implies that phenomena 


familiar in hydrodynamics should arise in gravity. In particular, motivated 
by this, arguments were presented that gravity could exhibit turbulent dy- 
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namics, with a direct energy cascade in 4 + 1 dimensions and the opposite 
in 3 + 1 |278j . Furthermore, in 3 + 1 dimensional gravity a quasi-conserved 
quantity should arise that is related to the conservation of entrophy in hy¬ 
drodynamics [279j . These observations have recently been demonstrated in 
ground-breaking numerical simulations of perturbed black branes |28n| (see 
Fig. 7.10), showing that the horizon geometry reflects the turbulent behav¬ 
ior of the boundary projection, and develops a fractal-like structure over the 
corresponding range of lengthscales. 


7 . 3 .7 Singularities 

Numerical simulations have played a significant role in analyzing the nature 
of singularities, in particular those which are often called “cosmological” due 
to the spacetimes having compact spatial topology (for a thorough review, 
see [283]). One of the longstanding questions has been the generic nature of 
singularities; i.e., what is the geometry of a spacetime approaching a singu¬ 
larity if no symmetries are presumed? For a vacuum spacetime, much of the 
research was inspired by the early work of Belinski, Lifschitz and Khalat- 
nikov (BKL,[283|), who conjectured that the generic singularity is spacelike, 
local, and oscillatory. The “local” part of the conjecture is that, in an appro¬ 
priate gauge, the spatial gradients in the field equations become irrelevant 
compared to the temporal gradients, and hence the dynamics at any spatial 
point reduces to a set of ordinary differential equations in time. The oscil¬ 
latory (or “mixmaster”) aspect then describes the dynamics of one of these 
points, claiming that the solution consists of an infinite, chaotic sequence of 
transitions between epochs, and in each epoch the geometry is well-described 
by one member of the Kasner family of geometries. A Kasner geometry is 
a homogeneous but anisotropic solution to the field equations consisting of 
two contracting and one expanding spatial direction (in the approach to the 
singularity). Several objections were raised to the BKL conjecture, in par¬ 
ticular that the assumptions they employed restricted their conclusions to 
local aspects of homogeneous cosmologies, and hence had little bearing on 
the generic, global properties of the spacetimes |285| . Numerical simulations 
have been key in resolving these disputes, gathering evidence in favour of 
the BKL conjecture [38t 128611287] . though discovering a surprising caveat in 
the process. This discovery was of so-called “spikes” that develop at isolated 
regions in the geometry [38] . A spike is a small lengthscale feature where the 
spatial gradients are not small, and hence are important in governing the lo¬ 
cal dynamics of the geometry. Spikes seem to undergo oscillatory transitions 
similar to the mixmaster behavior of non-spike worldlines |288j . However 
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Figure 7.9 Top: Energy density in planar shock collisions as a function of 
time t and longitudinal position z. The shocks approach each other along 
the z axis and collide at t = 0, z = 0. The collision produces “debris” that 
fills the forward light cone (from [7B].) Bottom: Depiction of the energy 
density of a 4d boundary flow dual to the evolution of a highly perturbed 
5d black hole in asymptotically global AdS spacetime (the radius of the 
black hole settles to 5 in geometric units, where the AdS length scale is 
L = \) (from |281| ). The boundary has topology and y is an an¬ 

gular coordinate; hence the image represents an initial high density (hence 
pressure) enhancement on the equator (y = 7r/2) that propagates back and 
forth between the equator and the poles (y = 0,7r). This result is from a 
pure 5d vacuum gravity simulation, yet the projected boundary dynamics 
matches that of a relativistic conformal fluid to within better than 1%, 
even in the early stages when the perturbation is highly non-linear. 
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Figure 7.10 Left: Vorticity of gravitational perturbations of a planar black 
hole as obtained through a 3 + 1 simulation of Einstein equations in AdS 
(from [76]). Right: Vorticity of a hydrodynamical field obtained in a 2 + 1 
viscous hydrodynamic simulation with a background fluid configuration 
dual to a planar black hole (from 1282) 1. Exploiting the fluid/gravity duality 
allows for constructing the full metric of the dual 3+1 spacetime to excellent 
accuracy. 

since (in the approach to the singularity) they shrink rapidly with time, 
even numerical simulations imposing planar symmetry (so 1 + 1 dimensional 
evolution) have not been able to follow their dynamics for long enough to 
conclusively demonstrate this. Due to these resolution challenges spikes have 
not been studied in scenarios with less symmetry, and so whether spike-like 
features beyond co-dimension 1 exist is also not known. 

An important point to make with regards to the above discussion of gener- 
icity and singularities is that it strictly applies only to these so-called cos¬ 
mological singularities, and not necessarily to those formed in gravitational 
collapse to black holes. There is some expectation that local properties of 
the singularities should be the same whether in a cosmological or black hole 
setting (indeed, the interior geometry of Schwarzschild is locally Kasner). 
However, the interior (Cauchy) horizon of rotating and charged black holes 
develop into null singularities when perturbed, and have very different struc¬ 
ture from the spacelike singularities in cosmology [289) . There are also argu¬ 
ments that null singularities are “as generic” as spacelike singularities [290) . 
and may also be relevant in a cosmological setting |291) . 


7.3.8 Miscellaneous 

Here we briefly discuss two miscellaneous topics where numerics have played 
an important role, and do not naturally fit into the main topic sections above. 
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Stability of AdS 

As opposed to Minkowski and deSitter spacetimes where global non-linear 
stability with respect to small perturbations has been established [29211293| , 
the related question has yet to be resolved in AdS. A key difference be¬ 
tween AdS and the other two spacetimes is the fact that infinity is timelike, 
and acts like a confining boundary; namely the future light cone from any 
event on an interior timelike observer’s worldline will reach the boundary 
and return to intersect the worldline a finite proper time later. Ground¬ 
breaking numerical and analytical work |ir)5| studied the spherically sym¬ 
metric Einstein-Klein-Gordon system in asymptotically AdS spacetime, and 
uncovered that a black hole eventually forms from an arbitrarily small ini¬ 
tial perturbation. This can heuristically be understood as a direct result of 
the confining property of AdS—energy cannot dissipate away, and due to 
non-linear interaction eventually a configuration will be explored where the 
central energy density becomes sufficiently large to cause gravitational col¬ 
lapse. A more quantitative explanation was given in uns], where through a 
resonant mechanism there is a secular transfer of energy from large to small 
scales, ending when a black hole forms. Furthermore, at the threshold of 
black hole formation, the spacetime behaves self-similarly, and the solution 
corresponds to the one seen in the asymptotically flat case [9l] (as expected 
since the AdS scale is irrelevant for a small black hole). Related work has 
argued that this behavior should still be present in the absence of symme¬ 
tries, and also when only gravitational perturbations are considered |294j . 
While these studies suggested AdS is unstable to arbitrarily small, generic 
perturbations, more recent follow up work has demonstrated the existence 
of large classes of initial data that are stable |29511296^ I297j . Applying these 
results to the AdS/CFT correspondence, given that black hole formation is 
synonymous with thermalization, this implies (perhaps unsurprisingly) that 
there are large classes of states in the dual GFT that do not thermalize. 


Formation and evaporation of CGHS Black Holes 

Two dimensional dilaton models of black hole evaporation were a popular 
subject of research a couple of decades ago, and though much was learned 
about the quantum nature of black holes from them, one could argue that 
no consensus results were obtained regarding the final near-Planck stages 
of evaporation (whether the black hole evaporates completely, if there is 
a remnant, a naked singularity, baby universe, etc.), or whether informa¬ 
tion is lost. One such popular model is that of Callan, Giddings, Harvey 
and Strominger (GGHS) [298]. Though extensively studied before, many 
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interesting quantitative and qualitative features of solutions to the semi- 
classical CGHS equations of motion were missed until a recent numerical 
study [2991130011301| . One of the more interesting results revealed here is 
that there are two distinct classes of solutions: those that can be identi¬ 
fied as microscopic (with initial masses of order the Planck mass or less), 
and those that are macroscopic (with initial masses a few times the Planck 
mass or larger). Remarkably, for macroscopic cases, after a brief transient, 
the evaporating spacetime and Hawking flux asymptote to a universal solu¬ 
tion, irrespective of details of the matter distribution that formed the black 
hole. Evaporation continues until the dynamical horizon shrinks to zero area, 
whence it encounters a singularity of the semi-classical equations (though 
this singularity is weaker than that arising in the classical solution). The 
future Cauchy horizon of this singularity is regular, in contrast with ear¬ 
lier suggestions that it would propagate to infinity in a “thunderbolt”. An 
improvement to the Bondi mass of the spacetime proposed in |302| shows 
that there is still on the order of a Planck mass “remnant” in the singular¬ 
ity, though this would presumably be resolved with higher order quantum 
corrections. 

This behavior is very different from that of black holes initially formed 
with only of the order of the Planck mass (the microscopic branch). Earlier 
studies had missed this distinction, and focused all attention on the physi¬ 
cally less relevant microscopic solutions. The macroscopic branch also turns 
out to be quite challenging to solve numerically, where scales of the order 
Mpianck in file initial vacuum state are exponentially “inflated” to scales 
of order g^/^pianck outgoing Hawking flux (a manifestation of the 

red-shift of outgoing radiation in black hole spacetimes, though in the evap¬ 
orating case this red-shift remains finite, since what was a null event horizon 
becomes instead a time-like dynamical horizon). 

One unusual aspect of 2d dilaton gravity that prevents straight-forward 
application of 2d results to the more relevant 4d case, is that in the former 
case there are two distinct, causally disconnected null infinities (“left” and 
“right”). This effectively disassociates the quantum state of the “ingoing” 
(right to left moving quanta, say) matter that forms the black hole from the 
“outgoing” (left to right) vacuum that becomes the Hawking flux. The semi- 
classical results together with the arguments presented in [302] suggest that 
the evolution of this vacuum sector is unitary, though little information 
about the infalling matter is retrievable from the Hawking flux. There is 
also no sign of any “firewall” [303) along the dynamical horizon at the semi- 
classical level. 
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Within the eternal inflation paradigm, our observable universe is contained 
in one of many bubbles formed from an inflating metastable vacuum |3n4| . 
Collisions between bubbles can potentially leave a detectable imprint on 
the cosmic microwave background radiation (see reviews |305t I306j ). While 
this scenario was initially studied through phenomenological models, recent 
works have concentrated on providing a quantitative connection between 
particular scalar field models giving rise to eternal inflation and the detailed 
signatures imprinted on the CMB. To this end, the intrinsically non-linear 
nature of the bubbles and their collisions have been studied numerically 
within full general relativity 130711308] , Simulations have revealed, in partic¬ 
ular, the following: i) the energy released in the collision of identical vacuum 
bubbles goes mostly into the formation of localized field configurations such 
as oscillons; ii) the structure of the potential considered is the dominant 
factor determining the immediate outcome of a collision; and iii) slow-roll 
inflation can occur to the future of a collision. Interestingly, these studies 
indicate that the signature in the CMB is well-described by a set of four 
phenomenological parameters whose values can be only probabilistically de¬ 
termined. 


Inhomogeneity in cosmology 

The majority of applications of general relativity to cosmology over the past 
decades have utilized analytical methods. For observational cosmology this 
is because the observed homogeneity and isotropy of the universe implies 
that its large scale structure is well-described by known exact solutions (the 
Friedmann-Robertson-Walker-Lemaitre metrics), with deviations from the 
FRWL solutions small and hence amenable to treatment by perturbation 
theory. There has however been some concern, in particular in light of the 
discovery of the present day accelerated expansion of the universe, that large 
scale inhomogeneities such as filaments and voids, or small scale non-linear 
inhomogeneities such as stars, can alter the assumptions made to study 
the largest scale dynamics of the cosmos (see [309] for a recent review). 
Some of these questions can be addressed by numerical solutions within full 
general relativity, in particular whether local, non-linear inhomogeneities can 
affect the overall expansion compared to a homogeneous universe with the 
same average stress-energy context (though “averaging” is itself an issue 
of some delicacy here). Only recently have simulations of such scenarios 
been considered |3in| 1311312] . In the latter two studies, universes with 
a positive cosmological constant and filled with a periodic lattice of black 
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holes (thus the most extreme example of non-linearity possible in general 
relativity) were evolved. The results were that the effective expansion rate 
was consistent with that of an equivalent homogeneous dust-filled universe. 


7.4 Unsolved Problems 

The aforementioned list of studies, while impressive in its own right, is only a 
portion of interesting phenomena where numerical relativity can shed light 
on important questions, as well as open up new research directions from 
them. The following is a (necessarily incomplete) list of such questions. 

• Strongly gravitating/highly dynamical scenarios and astrophysics. While 
it is clear that simulations have played a key role in uncovering the be¬ 
havior and characteristics of gravitational wave emission from compact 
binaries in astrophysical settings, much work, and many opportunities, 
remain. Indeed, even in the case of binary black holes where “only” the 
Einstein equations in vacuum are required, higher mass ratios and/or 
nearly-maximal spinning configurations have proven difficult and costly. 
The possible existence of intermediate mass black holes strongly motivates 
understanding the former class of binary. Non-vacuum systems require a 
more complex description due to the additional, and often involved, mat¬ 
ter physics. The rewards for this complexity are that now in addition 
to gravitational waves, electromagnetic and/or neutrino emission become 
possible, with the consequence that the opportunity for simulations to 
make contact with observation is extremely rich. The overarching goals of 
such simulations are to obtain first-principles descriptions of the detailed 
observational signatures across the range of emission channels the binaries 
might produce. The challenge to do so comes largely from the disparate 
time/length scales introduced by a plethora of physical processes, and 
by the complexity of the microphysics. Unlike the Einstein equations, to 
make simulations of realistic matter tractable invariably requires simpli¬ 
fied models of the fundamental equations. There is much opportunity here 
for synergy among the relevant communities: numerical relativity, gravi¬ 
tational wave observation, theoretical and observational astronomy, and 
nuclear physics. 

• Fundamental questions. Numerical simulations will continue to play a key 
role in exploring questions about the fundamental nature of Einstein grav¬ 
ity. There is no shortage of tantalizing questions remaining to be explored, 
including the non-linear development of superradiant and other black hole 
instabilities in four and higher dimensional spacetimes (see [313] for a first 
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such study in 4d), the nature of generic singularities inside rotating black 
holes (c.f. the “mass inflation” phenomenon |289j l. the dynamics of near¬ 
extremal black holes, cosmological domain wall and gravitational shock- 
wave collisions without symmetry assumptions, critical collapse without 
symmetries, black hole collisions and other dynamical non-linear inter¬ 
actions in asymptotically AdS spacetimes (in particular with AdS/CFT 
applications in mind), testing the limits of the hoop and cosmic censorship 
conjectures, the possible development -and consequences- of turbulence 
in gravity and fractal horizon structures, etc. (some of these are discussed 
further below). If past discoveries, such as critical phenomena and the 
“turbulent” instability in AdS spacetimes are any indication, many sur¬ 
prises await. Furthermore, these might have counterparts in other phys¬ 
ical systems, and important insights might be gained in both directions 
by analogy and their mathematical similarity or equivalence. 

Critical collapse. As reviewed above, the vast majority of studies of the 
threshold of gravitational collapse have been carried out in spherical sym¬ 
metry. That the original study of the axisymmetric gravitational wave crit¬ 
ical solution has defied attempts at a detailed solution for almost twenty 
years now hints at a very interesting and rich geometric structure await¬ 
ing discovery. For axisymmetric scalar field collapse, the inconsistency be¬ 
tween perturbative results suggesting that all non-spherical perturbations 
decay and a numerical study that hinted at a second, “focusing” instabil¬ 
ity remains to be resolved. Collapse without any symmetry assumptions 
is essentially uncharted territory. 

Black hole instabilities in higher dimensional, asymptotically flat space- 
times. A number of black holes in higher dimensions have been argued 
to be unstable, in particular by making connection to Gregory-Laflamme 
type instabilities. These arguments stem from the realization that the 
near-horizon geometry in (at least portions of) these black hole space- 
times can be mapped to unstable black string solutions, and so should 
display related phenomenology. This is the case for ultra-spinning black 
holes, black rings, black Saturns, etc. (e.g. [255]). Whether in all cases 
these black holes yield the rich behavior observed in perturbed, unstable 
black strings is yet unknown. For instance, in rapidly spinning black holes 
the instability induces a non-trivial, time-dependent quadrupole that ra¬ 
diates angular momentum that could shut off the instability. This has al¬ 
ready been observed in numerical simulations of Myers-Perry black holes 
in 6-8 dimensions, though only for cases with relatively mild angular mo¬ 
menta |258| . For sufficiently large angular momentum (recall that there 
is no upper bound in higher dimensions) the time scale of the Gregory- 
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Laflamme instability is shorter than the expected gravitational wave emis¬ 
sion time required to reduce the spin by enough to stabilize the system. A 
related problem is to consider highly prolate, “cigar-shaped” black holes 
in higher dimensions. Certainly, barring small scale length-wise perturba¬ 
tions, such a black hole would tend to becoming spherical on a time-scale 
of order equal to the light crossing time tct of the black hole. How¬ 
ever, such horizons that are sufficiently thin should locally be Gregory 
Laflamme unstable on a time scale much quicker than tct- Both of the 
aforementioned problems appear tractable in the near future. 

• Gravitational behavior in d ^ A and holography. As discussed above, 
holography has opened the door for numerical relativity to be exploited 
in problems outside the gravitational arena. Indeed, applications rele¬ 
vant to quark-gluon plasmas, condensed matter physics and quantum 
quenches (processes in which the physical couplings of a quantum sys¬ 
tem are abruptly changed), have recently been undertaken. While there 
is already an impressive body of work in this context, it is important to 
point out that most studies in this field have been in non-dynamical set¬ 
tings, and existing dynamical studies have assumed symmetries to yield a 
tractable computational problem. As a consequence, current results have 
certain limitations to the applicability and generality of the physics than 
can be drawn from them. This leaves much room for novel future work. 

• High speed/soliton collisions. Many questions remain in this topic. For 
soliton collisions, the nature of the black hole formation threshold solu¬ 
tion is unknown; possibilities include a “universal” gravitational critical 
solution irrespective of the nature of the matter, or alternatively the crit¬ 
ical solution of the matter field that the soliton is composed of. In the 
infinite boost limit, the geometry to the causal future of the shock wave 
collision is unknown. Very few studies of finite boost, higher dimensional 
collisions relevant to super-Planck scale particle collisions have been con¬ 
ducted. In particular only trivial topologies without brane tension have 
been considered, and charge has been ignored, which could be important 
at LHC energies. The intriguing suggestion of naked singularity formation 
in grazing 5d collisions shown in |229) needs further investigation. A de¬ 
tailed study of the radiation emitted in large impact parameter encounters 
in 4d would allow comparison with the effective field theory calculations 
that suggest the problem simplifies in this limit |239t I240j (the difficulty 
with such a study is that the black holes loose little energy in the en¬ 
counter, and hence a very long time numerical evolution will be required 
to allow the gravitational waves to get sufficiently far ahead of the black 
holes, unless novel gravitational wave extraction methods are developed). 
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For high speed collision applications to heavy ion collisions via AdS/CFT, 
future work includes relaxing symmetries to model non-central collisions, 
and introducing refinements to allow the dual CFT to better approximate 
QCD (for example, trying to model effects of confinement with additional 
matter fields, or via dynamics in the manifold of AdS^ x that are 
usually assumed to be trivial). 

• Alternative theories of gravity. Numerical relativity has also recently ven¬ 
tured into studying astrophysical binary systems within alternative grav¬ 
ity theories. Incipient investigations within Scalar-Tensor theories have 
uncovered an unexpected dynamical scalarization phenomena driven by 
the dynamics of binary neutron stars [314] . This phenomena has a signifi¬ 
cant impact on the orbiting behavior with clear consequences for the gravi¬ 
tational wave signals from the system. In all likelihood this behavior is only 
a token of the rich phenomenology awaiting to be discovered upon closer 
inspection of relevant theories, and can have astrophysical/observational 
consequences. 

• Supernova. Most core-collapse supernova simulations to date have not in¬ 
corporated full general relativity. Given that the problem of the explosion 
mechanism(s) is still unsolved and likely to be highly sensitive to the un¬ 
derlying physics, making the codes fully relativistic is another crucial step 
in the direction of more realistic modeling of the physics of 

this highly complex problem. 

As is clear from this list, there is no shortage of interesting applications 
for numerical studies. With the rapid development of numerical relativity 
over the past decades and its expansion to fields outside of pure classical 
general relativity, it is impossible to tell what a future review might have 
in store. At the same time however, it is safe to predict that many exciting 
results will fill its pages! 
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